Energetic particle tracing in optimized quasi-symmetric stellarator equilibria

Abstract

Recent developments in the design of magnetic confinement fusion devices have allowed the construction of exceptionally optimized stellarator configurations. The near-axis expansion in particular has been proven to enable the construction of magnetic configurations with good confinement properties while taking only a fraction of the usual computation time to generate optimized magnetic equilibria. However, not much is known about the overall features of fast-particle orbits computed in such analytical, yet simplified, equilibria when compared with those originating from accurate equilibrium solutions. This work aims to assess and demonstrate the potential of the near-axis expansion to provide accurate information on particle orbits and to compute loss fractions in moderate to high aspect ratios. The configurations used here are all scaled to fusion-relevant parameters and approximate quasi-symmetry to various degrees. This allows us to understand how deviations from quasi-symmetry affect particle orbits and what are their effects on the estimation of the loss fraction. Guiding-centre trajectories of fusion-born alpha particles are traced using gyronimo and SIMPLE codes under the NEAT framework, showing good numerical agreement. Discrepancies between near-axis and magnetohydrodynamic fields have minor effects on passing particles but significant effects on trapped particles, especially in quasi-helically symmetric magnetic fields. Effective expressions were found for estimating orbit widths and passing–trapped separatrix in quasi-symmetric near-axis fields. Loss fractions agree in the prompt losses regime but diverge afterwards.

1. Introduction

Energetic alpha particles generated by fusion reactions carry a significant amount of energy and have the potential to drive a plasma towards a self-sustaining fusion state, commonly referred to as a burning plasma. To achieve this, the fraction of alpha particles that deposit their energy in the plasma before being expelled needs to be maximized to guarantee sufficient alpha heating (Freidberg 2007). Additionally, energetic particles that leave the plasma can cause significant damage to the plasma-facing components of the fusion device, leading to a shorter device lifetime (Mau et al. 2008). Therefore, it is crucial to confine these energetic particles and accurately predict their behaviour in the plasma to achieve the desired levels of alpha heating and advance the development of fusion energy.

The most general condition for the confinement of orbits in stellarators is omnigenity, which requires a vanishing time-averaged radial drift, $\langle \Delta \psi \rangle$. An important subset of omnigenity, quasi-symmetry (QS), is obtained through a symmetry in the modulus of the magnetic field $\boldsymbol {B}$ when expressed in Boozer coordinates (Boozer 1981). Experiments such as the W7-X (Beidler et al. 1990) and HSX (Anderson et al. 1995) intend to approximate omnigenity via quasi-isodynamic and quasi-symmetric fields, respectively. Although perfectly omnigenous (Cary & Shasharina 1997) and quasi-symmetric (Garren & Boozer 1991a) devices have been conjectured not to exist, it has been shown that very precise approximations can be obtained (Landreman & Paul 2022; Goodman et al. 2023).

While these experiments have been proven to generally confine thermal plasmas, the loss of energetic alpha particles is still a key research area for stellarators (Bader et al. 2019a; LeViness et al. 2022). Minimizing the amount of lost particles is usually performed through optimization of plasma properties such as quasi-symmetry (Henneberg et al. 2019) and omnigenity (Goodman et al. 2023) or using proxies for loss fractions, such as $\varGamma _c$ (Bader et al. 2019b) and $\varGamma _{\alpha }$ (Sánchez et al. 2023). Direct optimization of loss fractions has only recently been performed, albeit at a substantial computational cost (Bindel, Landreman & Padidar 2023).

Such optimization procedures require the repeated variation of the last closed flux surface and are not only hindered by the existence of multiple local minima, and thus highly dependent on the initial conditions, but are also computationally demanding and offer limited insight into the number of effective degrees of freedom of the problem, as well as the specific physical implications associated with individual coefficients or their collective groups (Landreman, Sengupta & Plunk 2019). Additionally, to obtain a reactor-grade stellarator, it is essential to balance the confinement of energetic particles with other physical and engineering constraints such as stability, turbulence and coil geometry tolerances (Hegna et al. 2022), which may increase its overall computational cost.

A way to circumvent the issues above is through the introduction of an analytical approximation to a magnetohydrodynamic (MHD) equilibrium. Such construction is a near-axis expansion (NAE), which is valid in the core of all stellarators and can be introduced in the early stages of optimization enabling the use of better initial conditions (warm start initial conditions) for conventional optimizations while allowing for extensive searches in the parameter space of the design (Landreman 2022). This expansion can be obtained resorting either to the direct method (Mercier 1964; Solov'ev & Shafranov 1970; Jorge, Sengupta & Landreman 2020), explicitly determining the magnetic flux surface function $\psi$ using the Mercier coordinates ($\rho$, $\theta$, $\varphi$), or the indirect method (Garren & Boozer 1991a,b; Landreman & Sengupta 2018; Landreman et al. 2019), which directly computes the position vector $\boldsymbol {r}$ as a function of the Boozer coordinates and is the one we will use throughout this work. Such constructions not only describe high aspect ratio devices but also the region around the axis of any stellarator, while enabling useful analytical insight due to its simplicity at lowest orders. As we show here, the estimation of loss fractions directly from orbit following codes can be computationally expedited within the near-axis framework, without significant accuracy degradation due to its simplified analytical nature for short time frames. This speed up arises from the direct computation of the magnetic field and other related quantities as an alternative to the many Fourier coefficients needed in conventional optimization.

The main goal of this paper is to study the accuracy of fast particle trajectories in approximate analytical near-axis equilibria when compared with the same trajectories in accurate numerical MHD equilibria. By understanding how much these orbits deviate from each other, we are able to assess the validity of the expansion for loss fractions estimation in scenarios with a variety of aspect ratios, time scales and perpendicular velocities. As the collisional effects have been shown to be negligible for particle losses in the prompt losses time scales and time scales up to a fraction of the energy deposition time (Lazerson, LeViness & Lion 2021), we are only interested in the collisionless dynamics of the energetic particles, which were studied resorting to the integration of the particles’ guiding-centre (GC) orbits. The primary focus of this work is on configurations that approximate quasi-symmetry to varying degrees to employ the quasi-symmetric version of the NAE of Landreman & Sengupta (2019) and Landreman (2022). The concept of quasi-symmetry, derived from the gyroaveraged Lagrangian, implies that QS stellarators are specifically designed to confine GC orbits. Therefore, our analysis focused solely on these orbits. This approach is also more computationally efficient than simulating the full orbits of particles.

A large number of GC or full-gyromotion particle tracers, such as ANTS (Drevlak et al. 2014), ASCOT (Hirvijoki et al. 2014), BEAMS3D (McMillan & Lazerson 2014), FOCUS (Clauser, Farengo & Ferrari 2019), LOCUST (Ward et al. 2021), OFMC (Tani et al. 1981), GNET (Masaoka & Murakami 2013), SPIRAL (Kramer et al. 2013) and VENUS-LEVIS (Pfefferlé et al. 2014), have been used for fast particle transport and loss studies. In this work, we use the Euler–Lagrange equations of motion for the GC (Littlejohn 1983) as implemented in the general-geometry particle tracing library gyronimo (Rodrigues et al. 2023), which is a convenient tool for comparing between different geometries, and the Hamilton equations of motion as implemented in the code SIMPLE (Albert, Kasilov & Kernbichler 2020b). Both particle-tracing codes are open source. To make the direct comparison between trajectories as straightforward as possible, the same equations and algorithms are used for particles in near-axis and MHD fields. Later specialization could improve the computational efficiency of the procedures. However, the achieved levels of speed are already relevant for near-axis equilibria optimization.

This paper is organized as follows. Section 2 provides an overview of the NAE and the Boozer coordinate system. In § 3, we describe the GC motion and the different approaches taken for its calculations. In § 4, the results for single-particle tracing in configurations close to quasi-symmetry are shown, within a reasonable scope of initial conditions. Additionally, expressions for the estimation of the orbit radial amplitudes and the passing–trapped boundary are derived and compared with computational results. Taking knowledge obtained in this last section into account, total loss fractions are computed and analysed in § 5 for two distinct configurations. In § 6, we summarize the primary findings and delineate potential pathways for future research.

2. The near-axis expansion

In this section, we follow the notation of Landreman & Sengupta (2019) and introduce the near-axis expansion using the inverse approach first presented by Garren & Boozer (1991a). We begin by writing the magnetic field in Boozer coordinates $(\psi,\theta,\varphi )$, with $2{\rm \pi} \psi$ the toroidal flux, and $\theta$ and $\varphi$ the poloidal and toroidal coordinates, respectively. This results in

(2.1)\begin{align} \boldsymbol{B} & =\boldsymbol{\nabla}\psi \times\boldsymbol{\nabla}\theta + \iota \boldsymbol{\nabla}\varphi \times\boldsymbol{\nabla}\psi, \nonumber\\ & =\beta \boldsymbol{\nabla}\psi + I \boldsymbol{\nabla}\theta + G \boldsymbol{\nabla}\varphi, \end{align}

where $I=I(\psi )$ and $G=G(\psi )$ are flux functions, $\beta =\beta (\psi,\varphi, \theta )$ is a quantity related to the plasma pressure that usually depends on the three coordinates and $\iota =\iota (\psi )$ is the rotational transform. As we want to analyse both stellarators which exhibit a quasi-symmetry in the toroidal angle, named quasi-axisymmetric (QA), and in a linear combination of poloidal and toroidal angles, termed quasi-helically symmetric (QH), it is convenient to introduce a helical angle $\vartheta = \theta - N \varphi$, where $N$ is a constant integer equal to 0 for a QA stellarator and, usually, the number of field periods for a QH one. This results in the following expression for the magnetic field:

(2.2)\begin{align} \boldsymbol{B} & = \boldsymbol{\nabla}\psi \times\boldsymbol{\nabla}\vartheta + \iota_N \boldsymbol{\nabla}\varphi \times\boldsymbol{\nabla}\psi, \nonumber\\ & =\beta \boldsymbol{\nabla}\psi + I \boldsymbol{\nabla}\vartheta + (G+NI) \boldsymbol{\nabla}\varphi. \end{align}

To derive the NAE's main features at first order, we begin by writing the general position vector $\boldsymbol {r}$ as

(2.3)\begin{equation} \boldsymbol{r}(r,\theta, \varphi) = \boldsymbol{r}_0(\varphi) + r X_1(\theta, \varphi ) \boldsymbol{n}(\varphi) + r Y_1(\theta, \varphi ) \boldsymbol{b}(\varphi) + \mathcal{O}(r^2), \end{equation}

with $r = \sqrt {2|\psi |/\bar {B}}$ an effective minor radius, $\bar {B}$ a reference magnetic field and $\boldsymbol {r}_0$ the position vector along the axis, with a given set of Frenet–Serret orthonormal basis vectors ($\boldsymbol {t}$, $\boldsymbol {n}$, $\boldsymbol {b}$), a local curvature $\kappa =\kappa (\varphi )$ and torsion $\tau =\tau (\varphi )$. At first order, we additionally take the profile functions $G(r)$ and $I(r)$ to be $G_0$ and $I_2 r^2$, respectively (Landreman & Sengupta 2018). In quasi-symmetry, $\boldsymbol {r}$ can be written as (Landreman & Sengupta 2019)

(2.4)\begin{align} \boldsymbol{r}(r,\vartheta,\varphi) = \boldsymbol{r_0}(\varphi) + \frac{r \bar{\eta}}{\kappa(\varphi)}\cos\vartheta \boldsymbol{n}(\varphi) + \frac{r s_\psi s_G \kappa(\varphi)}{\bar{\eta}} [\sin\vartheta + \sigma(\varphi) \cos\vartheta] \boldsymbol{b}(\varphi) +O(r^2/\mathcal{R}), \end{align}

where $s_\psi = \mathrm {sign}(\psi )$, $s_G = \mathrm {sign}(G_0)$ and $\bar {\eta }$ is a constant reference field strength that parametrizes $B= |\boldsymbol {B}|$ as

(2.5)\begin{equation} B = B_0 (1 + r \bar{\eta} \cos\vartheta) + O((r/\mathcal{R})^2), \end{equation}

where $\bar {B}=s_\psi B_0$ and $\sigma (\varphi )$ as a periodic function related to the flux surface shape, which satisfies the Riccati-type equation

(2.6)\begin{equation} \frac{{\rm d}\sigma}{{\rm d}\varphi} + (\iota_0 - N) \left[ \frac{\bar{\eta}^4}{\kappa^4} + 1 + \sigma^2 \right] -\frac{2 G_0 \bar{\eta}^2}{B_0 \kappa^2} \left[ \frac{I_2}{B_0} - s_\psi \tau \right] = 0. \end{equation}

Equation (2.5) only enables the construction of stellarators with elliptical cross-sections, so a higher order is needed to express the stronger shaping of existing stellarators. At the second order, nine new functions of $\varphi$ arise in the surface shapes that can now possess triangularity and a Shafranov shift. However, these functions are constrained by 10 new equations of $\varphi$, a mismatch that results in the fact that most axis shapes are not consistent with quasi-symmetry at this order. Landreman & Sengupta (2019) circumvented this by allowing quasi-symmetry to be broken at second order in the field strength, which is the method used here. Furthermore, a detailed description of the application of this method to the construction of stellarator shapes is provided, built on a third-order method with the shaping details that need to be taken into account to generate a boundary surface that is consistent with the desired field strength and eliminate unwanted mirror modes.

As the present work performs comparisons between equilibria generated with the quasi-symmetric NAE code pyQSC (Landreman 2021) and with VMEC (Hirshman 1983), an MHD equilibrium code, it is important to note that these do not use the same coordinates. Unlike the description made for the radial coordinate in the NAE, VMEC uses a normalized radial coordinate $s=\psi /\psi _b$, with $\psi _b$ the toroidal flux at the last closed flux surface. Therefore, to obtain $r$ in the NAE system, we use the relation

(2.7)\begin{equation} r = r_{\text{max}} \sqrt{s}, \end{equation}

where $r_{\text {max}}=\sqrt {2 \psi _b /\bar {B}}$. Additionally, the toroidal angle in VMEC, $\phi$, corresponds to the azimuthal cylindrical coordinate defined as $\arctan (\kern0.7pt y/x)$, while pyQSC operates with a cylindrical coordinate on-axis, $\phi _0$, that corresponds to the asymptotic value on axis, i.e.

(2.8)\begin{equation} \phi_0= \lim_{r \to 0} \tan^{{-}1} \left(\frac{y}{x}\right),\end{equation}

which is distinct from the Boozer coordinate described in the NAE derivation above. In our workflow, we use the Boozer coordinate $\varphi$ and transform it into $\phi _0$ through $\phi _0 = \varphi - \nu$, where $\nu$ is an output of the pyQSC code given by Landreman & Sengupta (2018). The coordinate $\phi$ is further computed through a root-finding function. Finally, the poloidal coordinate in VMEC, $\theta _V$, is in the range $[-{\rm \pi},{\rm \pi} ]$, whereas $\theta$ in pyQSC is defined from $0$ to $2 {\rm \pi}$ in the opposite direction, leading to the relation $\theta _V={\rm \pi} - \theta$. The initialization of a particle in a pyQSC field is, however, done with the helical angle $\vartheta =\theta - N \varphi$.

3. Guiding-centre formalism and orbit integration

The confinement properties of a given stellarator configuration are commonly estimated employing either proxies or direct methods, such as numerically simulating particle orbits in a plasma. Although proxies for confinement such as QS, $\varGamma _c$ and others have proven fruitful in device optimization due to their minimal computational expense (Bader et al. 2019a; Henneberg et al. 2019), they also bring some drawbacks. For example, QS may be too stringent of a condition when dealing with a multi-objective optimization, possibly undermining our capacity to achieve optimal outcomes for other objectives like MHD stability or feasible coil sets, and proxies like $\varGamma _c$ and $\varGamma _\alpha$ might not capture the full physical picture of drift orbits (Albert, Kasilov & Kernbichler 2020a). Furthermore, orbit following is necessary if one wants to either understand the mechanisms that lead to losses or take those into account in the loss estimation. Thus, we apply the equations of motion for charged particles in a magnetized plasma to assess such mechanisms.

Computing particle trajectories in a magnetic field is commonly not a straightforward task, often requiring numerical treatment. To do so, we leverage the fact that in the magnetized plasmas of interest, the magnetic field is strong enough that we can decouple the high-velocity cyclotronic motion of charged particles from the slowly varying GC coordinates and parallel velocity, therefore simplifying our problem. We begin by writing the position of a particle $\boldsymbol {r}$ in terms of its GC position $\boldsymbol {R}$ and its gyration radius vector $\boldsymbol {\rho }$ in the following way (Cary & Brizard 2009):

(3.1)\begin{equation} \boldsymbol{r}(t) = \boldsymbol{R}(t) + \boldsymbol{\rho}(t). \end{equation}

If the gyroradius $\rho _c=mv_\perp /qB$ of a particle of mass $m$ and charge $q$ is small enough compared with a field gradient length $L_B$, we can average out the particle's motion over a gyration time $t_c=1/\varOmega _c$, where $\varOmega _c=qB/m$ is the gyrofrequency, and still retain most of the information about the particle motion for time scales much greater than $t_c$. To a first approximation, the GC motion can be described with the non-canonical gyro-averaged Lagrangian (Littlejohn 1983):

(3.2)\begin{equation} \mathcal{L} = q \boldsymbol{A} \boldsymbol{\cdot} \dot{\boldsymbol{R}} + m v_{{\parallel}}\boldsymbol{b} \boldsymbol{\cdot} \dot{\boldsymbol{R}} + \frac{m \mu}{q} \varOmega_c - \frac{m v_{{\parallel}}^2}{2} - \mu B, \end{equation}

where $\mu = m v_\perp ^2 / 2 B$ is an adiabatic invariant (magnetic moment), $v_\parallel$ and $v_\perp$ are respectively the components of the velocity parallel and perpendicular to the magnetic flux surfaces, $\boldsymbol {b} = \boldsymbol {B}/B$ the magnetic unit vector, and $\boldsymbol {A}$ is the magnetic potential. As we are not considering the effects of electric fields from local charge densities or varying magnetic fields, the electric potential is absent from the Lagrangian above. The application of Lagrangian mechanics to this system allows us to use the Euler–Lagrange equations to derive equations of motion that not only conserve energy for time-independent systems but also possess Poincaré integral invariants and allow for derivations of conservation laws from Nöether's theorem in the presence of spatial symmetries (Cary & Brizard 2009). The GC Lagrangian in Boozer coordinates has the following form:

(3.3)\begin{equation} \mathcal{L} = q (\psi \dot{\theta} -\psi_p \dot{\varphi}) + \frac{m v_{{\parallel}}}{B} (\dot{\psi} \beta + \dot{\theta} I + \dot{\varphi} G) + \frac{m \mu}{q} \varOmega_c - \frac{m v_{{\parallel}}^2}{2} - \mu B, \end{equation}

where we used the fact that, as $\boldsymbol {B}=\boldsymbol {\nabla } \times \boldsymbol {A} = \boldsymbol {\nabla } \psi \times \boldsymbol {\nabla } \theta + \boldsymbol {\nabla } \varphi \times \boldsymbol {\nabla } \psi _p$ with ${\rm d} \psi _p / {\rm d} \psi =\iota$, $\boldsymbol {A}$ can be written as $\boldsymbol {A}=\psi \boldsymbol {\nabla } \theta - \psi _p \boldsymbol {\nabla } \varphi$, where $\psi _p$ is the poloidal flux. In the case of quasi-symmetry, where the magnitude of $\boldsymbol {B}$ depends only on $(\psi,\chi )$, with $\chi =M \theta - N \varphi$, (3.3) exhibits a continuous symmetry in a third coordinate $\eta =M' \theta - N' \varphi$, with $M'/N' \neq M/N$. In compliance with Nöether's theorem, this symmetry leads to the conservation of a quantity analogous to the canonical angular momentum in the tokamak, which constrains the trajectories of the particles in the plasma, improving their overall confinement up to modern tokamak levels or greater (Landreman & Paul 2022).

Once we have the equations of motion, there is no unique method of performing their integration to obtain particle orbits. Although (3.2) is valid independently of the applied system of coordinates, the standard approach to particle tracing is to write this Lagrangian in terms of coordinates that are more suitable to the problem in question, such as Boozer coordinates, as shown in (3.3), usually leading to simplified equations of motion. This is the case of SIMPLE (Albert et al. 2020b), a symplectic (energy conserving) orbit tracer used in this work. The symplectic approach to particle following has significant advantages regarding numerical stability and energy conservation, but its implementation is not straightforward. Symplectic integrators rely on canonical coordinates and an explicit expression of the GC Hamiltonian,

(3.4)\begin{equation} H(\boldsymbol{z}) = \frac{m v_\parallel^2(\boldsymbol{z})}{2} + \mu B(\boldsymbol{z}), \end{equation}

which exists only in non-canonical coordinates $\boldsymbol {z}=(\psi,\theta,\varphi, p_\varphi )$, where $p_\varphi = \partial _{\dot {\varphi }} \mathcal {L}$, so a coordinate transformation to canonical coordinates $(\boldsymbol {q},\boldsymbol {p})=(\theta,\varphi, p_\theta,p_\varphi )$ is needed and can be found from Albert, Kasilov & Kernbichler (2020c), with $p_\theta = \partial _{\dot {\theta }} \mathcal {L}$. The canonical coordinates can be written as explicit and invertible functions $\boldsymbol {q} = \boldsymbol {q} (z)$ and $\boldsymbol {p} = \boldsymbol {p} (z)$ with inverses $\boldsymbol {z} (\boldsymbol {q}, \boldsymbol {p})$ given only implicitly. Although this would not be sufficient for an explicit integrator, it is enough for semi- or fully implicit integrators.

The equations of motion of the canonical coordinates in SIMPLE are given by the following set of equations:

(3.5)\begin{gather} \dot{\theta}(t) =\frac{\partial H}{\partial \psi}\left(\frac{\partial p_{\theta}}{\partial\psi}\right)^{{-}1}, \end{gather}

(3.6)\begin{gather}\dot{\varphi}(t) =\frac{\partial H}{\partial p_{\varphi}}-\frac{\partial H}{\partial \psi} \left(\frac{\partial p_{\theta}}{\partial \psi}\right)^{{-}1} \frac{\partial p_{\theta}}{\partial p_{\varphi}}, \end{gather}

(3.7)\begin{gather}\dot{p}_{\theta}(t) ={-}\frac{\partial H}{\partial\theta}+\frac{\partial H}{\partial \psi} \left(\frac{\partial p_{\theta}}{\partial \psi}\right)^{{-}1}\frac{\partial p_{\theta}}{\partial\theta}, \end{gather}

(3.8)\begin{gather}\dot{p}_{\varphi}(t) ={-}\frac{\partial H}{\partial\varphi}+\frac{\partial H}{\partial \psi}\left(\frac{\partial p_{\theta}}{\partial \psi}\right)^{{-}1}\frac{\partial p_{\theta}}{\partial\varphi}, \end{gather}

which are then solved with symplectic schemes, such as implicit Euler schemes, Verlet and implicit midpoint, estimating the time evolution of the radial coordinate $\psi$ in the process, as further described by Albert et al. (2020a).

The inconveniences related to canonical transformations and their inversions can be avoided by employing the Euler–Lagrange equations of motion that result from the Lagrangian in (3.2), as these are invariant regarding any specific choice of coordinates. This is the approach followed by one of the equations-of-motion models implemented in the library gyronimo, with the corresponding dynamical system being written in terms of the general curvilinear coordinates $\boldsymbol {X}$ as

(3.9)\begin{gather} \varOmega^{*} \frac{{\rm d}\boldsymbol{X}}{{\rm d}\tau} = \tilde{\varOmega} \tilde{v}_\parallel \boldsymbol{b} + \tilde{v}_\parallel^2 \tilde{\boldsymbol{\nabla}} \times \boldsymbol{b} - \boldsymbol{b}\times\left(\tilde{\boldsymbol{E}} - \tilde{v}_\parallel\partial_\tau\boldsymbol{b} - \frac{\tilde{\mu}\tilde{\boldsymbol{\nabla}}\tilde{B}}{2}\right), \end{gather}

(3.10)\begin{gather}\varOmega^{*} \frac{{\rm d}\tilde{v}_\parallel}{{\rm d}\tau} = (\tilde{\varOmega} \boldsymbol{b} + \tilde{v}_\parallel \tilde{\boldsymbol{\nabla}}\times\boldsymbol{b}) \boldsymbol{\cdot} \left(\tilde{\boldsymbol{E}} -\tilde{v}_\parallel\partial_\tau\boldsymbol{b} - \frac{\tilde{\mu}\tilde{\boldsymbol{\nabla}}\tilde{B}}{2}\right), \end{gather}

where $\varOmega ^{*} = \tilde {\varOmega } + \tilde {v}_\parallel (\boldsymbol {b} \boldsymbol {\cdot } \tilde {\boldsymbol {\nabla }} \times \boldsymbol {b} )$, the position $\boldsymbol {X}$ of the GC is normalized to a reference length $L_\text {ref}$, the time $\tau$ to $T_\text {ref}$ and the parallel velocity to $V_\text {ref}=L_\text {ref}/T_\text {ref}$. In these expressions, $\tilde {B} = B/B_\text {ref}$ is a normalized magnetic field magnitude, $\tilde {\boldsymbol {E}} = \tilde {E}_\text {ref} (\boldsymbol {E}/E_\text {ref})$ a normalized electric field and $\tilde {\varOmega } = \varOmega _\text {ref} \tilde {B}$ a normalized gyrofrequency. Moreover, $\tilde {\boldsymbol {\nabla }} = L_\text {ref} \boldsymbol {\nabla }$, $\tilde {E}_\text {ref} = \tilde {\varOmega }_\text {ref} (E_\text {ref} V_\text {ref}^{-1} B_\text {ref}^{-1})$ and the magnetic moment $\mu$ is normalized to the ratio $U_\text {ref}/B_\text {ref}$ and $U_\text {ref}$ is the kinetic energy corresponding to $V_\text {ref}$.

This enables the use of different kinds of magnetic equilibria in any given geometry. Using these equations of motion, the integration of the particles’ dynamics can be performed independently of the specific geometry and coordinate system of the problem, taking only the tensor metric $\boldsymbol {g}$ and either the co- or contravariant components of $\boldsymbol {B}$ as inputs. All magnetic configurations analysed in this work are equilibrium vacuum configurations with no electric fields resulting in $\partial _\tau \boldsymbol {b}=0$, $\boldsymbol {b} \boldsymbol {\cdot } \tilde {\boldsymbol {\nabla }} \times \boldsymbol {b} = 0$ and $\tilde {\boldsymbol {E}}=0$, effectively reducing the equations of motion to

(3.11)\begin{equation} \frac{{\rm d}\boldsymbol{X}}{{\rm d}\tau} = \tilde{v}_\parallel \boldsymbol{b} + \frac{\tilde{\mu}}{2 \tilde{\varOmega}} \boldsymbol{b}\times \tilde{\boldsymbol{\nabla}}\tilde{B} \end{equation}

and

(3.12)\begin{equation} \frac{{\rm d}\tilde{v}_\parallel}{{\rm d}\tau}={-} \frac{\tilde{\mu}}{2} \boldsymbol{b} \boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}}\tilde{B}. \end{equation}

The equations of motion (3.9) and (3.10) can be integrated with a wide assortment of algorithms made available by the ordinary differential equation (ODE) library boost::odeint (Ahnert & Mulansky 2011), either with constant or adaptive time step, which can be easily selected from within the gyronimo framework and allow the flexibility to choose the algorithm best adapted to the problem being solved. Henceforth, we will work with the Dormand–Prince (RK5) method from the Runge–Kutta family of ODE solvers. This decision is based in figure 1, where panel (a) shows the trade-off between the energy loss $\Delta E= (E-E_0)/E_0$, which is a proxy for the computational error since the energy should be conserved, and the computation time required by different constant time step algorithms to integrate up to the same final instant. In panel (b), one compares the best-performing constant-step algorithms with their adaptive-step counterparts when available.

Figure 1. Evaluation of the performance of different integration algorithms. (a) Constant-step Runge–Kutta (RK), Adams–Bashforth (AB) and Bulirsch–Stoer (BS) integrators. (b) Best performing constant-stepintegrators and their adaptive-step counterparts available in the boost::odeint library.

Although the Bulirsch–Stoer (BS) method and the fourth-order Runge–Kutta (RK4) constant-step methods exhibit better performances, they entail increased computational costs, which is a significant bottleneck when evaluating loss fractions for a wide range of magnetic configurations in optimization workflows. Even though the adaptive-step Runge–Kutta–Fehlberg (RK78 adapt) appears to accomplish this, the parameters that control the energy loss do not show a linear dependence and so were left to later statistical and systematic analysis. It is also interesting to note that there is a limit on how small we can make the step size and still decrease the relative energy loss, a limit that may be due to the numerical precision used in the code.

4. Single-particle motion

We begin by performing integration of the GC trajectories of alpha particles at an energy of $3.52$ MeV. All of the equilibria were scaled to have an effective minor radius $a_A=1.7044$ m, equal to that of the ARIES-CS reactor (Najmabadi et al. 2008), and major radius larger than ARIES-CS, $R_A=7.7495$ m, with a magnetic field of $5.3267$ T. Here, ARIES-CS is used as a baseline scale as it was one of the latest viability studies for a power plant-grade compact stellarator. However, as it is a compact device, its aspect ratio $A$, the ratio between its major radius $R_A$ and its minor radius $a_A$, may stretch the limit of applicability of the NAE to be effective all the way up to the boundary, so we take higher values of the major radius to obtain a better behaviour of this approximation, albeit at the cost of a less compact reactor. The effective minor radius in this work is given as an average of this quantity throughout the last closed surface, which can be written as ${\rm \pi} a_{A}^2 = 1/2 {\rm \pi}\int _0^{2{\rm \pi} } S(\phi )\,{\rm d}\phi$, where $S(\phi )$ is the area of the surface cross-section at each point of the cylindrical angle $\phi$ (Jorge & Landreman 2021). We note that the minor radius value is substantially larger than the alpha particles’ gyroradius as expected of this kind of fusion reactor.

Finally, we initialize the particles velocities with two parameters: $s_{v_\parallel }$, the sign of the initial parallel velocity, which can take the values $+1$ and $-1$, and $\lambda$, the ratio between the initial perpendicular kinetic energy $E_{\perp _i}= m v_\perp ^2 / 2$ and the total kinetic energy $E=m v^2 / 2$. The normalized adiabatic invariant is given by $\tilde {\mu }=\lambda (E/E_\text {ref}) (B_\text {ref}/B(r_i))$, where $r_i$ represents the initial position of the particle.

To benchmark gyronimo with SIMPLE, the same alpha particles were followed with both tracers for $1$ ms in the precise QA VMEC equilibrium of Landreman & Paul (2022) scaled to the conditions described above, except for the major radius, which is $R = 6 a_A$ to preserve the original aspect ratio of $A = 6$. Figure 2 exhibits the orbits of particle with $s_{v_\parallel }=1$ and initial position $(s,\theta, \phi )=(0.25, 0.1, 0.1)$. The angular coordinates $\theta$ and $\phi$ were chosen to be finite values to ensure these are not a special case where the orbits coincide, which could happen for $\theta =\phi =0$. Although this benchmark included the variation of all initial values, we only exhibit here a change in the value of $\lambda$ which was observed to be one of the most important factors in the alignment between integrators. In the presented results, the values of $\lambda$ are $0.95$, $0.97$ and $0.99$.

Figure 2. Comparison between gyronimo (vmec) and SIMPLE (simple) tracers for original precise QA, scaled to have ARIES-CS minor radius $a_A$ and $B_0$, the field at the plasma axis. With an initial position $(s,\theta, \phi )=(0.25, 0.1, 0.1)$ and $\lambda =0.95$ in panels (a,d), the particle describes a passing orbit, and for $\lambda =0.97$ in panels (b,e) and $\lambda =0.99$ in panels (c, f), trapped orbits are observed. Top panels are the temporal evolution of the radial coordinate $s$ and bottom panels are the poloidal view of the orbit in VMEC coordinates, with an inner dashed circle for the initial flux surface and an outer dashed circle for the last closed flux surface. Only approaching the passing–trapped separatrix ($\lambda =0.97$) do the tracers slightly diverge.

The radial oscillations predicted by the two different approaches agree for the cases $\lambda =0.95$ and $\lambda =0.99$, which correspond respectively to a passing and a trapped particle. Such agreement, however, starts to decrease over extended temporal intervals for the case $\lambda =0.97$. This corresponds to a trapped particle that starts to approach the passing–trapped transition where orbits tend to be more sensitive to slight changes in the integration conditions. For those values of $\lambda$, the amplitude of oscillations in $s$ remains the same but its temporal variation contains a delay between the codes, while the average radial position exhibits a small deviation. While the temporal deviation may be due to the symplectic nature of SIMPLE, the source of the difference in the average radial position, among other factors, may include the fact that the motion of a particle in the transition from the trapped to the passing regime depends strongly on higher order terms of the GC motion. The segment of the banana orbits around reflection points (banana tips) correspond to low values of $v_\parallel$. We note that such higher order terms may have a non-negligible effect in points along the orbit where the values of $v_\parallel$ stay close to zero.

We conclude that we may work with gyronimo as it compares well in most cases with SIMPLE. Throughout the rest of the analysis, we will not be showing the results of the SIMPLE particle tracing as such a good match for some orbits is only achievable using an increased amount of memory due to the number of grid points necessary for the interpolation used by SIMPLE. The effect of varying an interpolation angular grid factor that controls the number of poloidal and toroidal grid points within SIMPLE can be seen in figure 3 for a particle initialized with $s_i=(0.25, 2.89, 1.84)$ and $v_\parallel$ $/ v = 0.44$, with drastic changes on its passing orbit. This is partially explained by the fact that the main goal of SIMPLE is to achieve long-term conservation of invariants like energy rather than minimizing the spatial deviation from a reference orbit (Albert et al. 2020c). Despite that, all single-particle results were compared with the corresponding values in SIMPLE with a low grid factor to ensure that there are no major divergences between tracers, including those shown in this work.

Figure 3. Orbits obtained with SIMPLE with $(s,\theta, \phi )=(0.25, 2.89, 1.84)$ and $v_\parallel / v = 0.44$ as initial conditions, by varying the angular grid factor multharm: (a) multharm $=3$; (b) multharm $=4$; (c) multharm $=5$.

We will divide our comparison for QA and QH configurations due to the nuanced analysis needed between them. The reference case for a QA stellarator used was the precise QA of Landreman & Paul (2022) and the one for QH stellarators was the four field period QH vacuum configuration with magnetic well from § 5.4 of Landreman (2022). Although we will use the terms QH and QA for all of the magnetic configurations that are close to precise QS, some have finite deviations from it. This is the case because we take the outer surface of NAE configurations calculated by pyQSC and create an input file for VMEC with the relevant configuration parameters, following the procedure of Landreman et al. (2019). VMEC computes an ideal MHD equilibrium from this input file which may not retain some of the initial features of the NAE configuration. This is especially true for equilibria with smaller aspect ratios, where the outer bounds are not well described by an expansion from the axis.

To have a good overall picture of how the NAE compares with with different aspect ratio equilibria, the analysis will encompass equilibria with major radius $R=3, 2$ and $1.5$ times the major radius of ARIES-CS, $R_A$. We can see from figures 4 and 7 how quasi-symmetry degrades when generating a VMEC equilibrium from the near-axis one for lower aspect ratios. While this effect is more evident for the QH stellarator in figure 7, such degradation is visible also in the QA one in the contour lines representing higher values of $B$ in figure 4. The contour lines representing lower magnitudes of $B$ appear not to follow this trend, but only because they are closer to $\theta / 2 {\rm \pi}= 0.5$ for the $A=13.6$ equilibrium than for the $A=9.1$ one.

Figure 4. QA VMEC equilibria generated from the precise QA from pyQSC with different major radius scalings. Top panels are the 3-D versions of the equilibria and bottom panels are the contour plots of the magnetic fields on the angular VMEC coordinates. (a) Aspect ratio $A=6.8$. (b) $A=9.1$. (c) $A=13.6$.

4.1. Quasi-axisymmetry

We now proceed to the comparison between orbits of particles for the VMEC equilibria and the pyQSC ones that originated them. All orbits were traced for $1$ ms, to be compared in the same time frame, although varying the aspect ratio affects the number of toroidal turns and the number of radial oscillations, similarly for changes in the value of $\lambda$.

We show in figures 5 and 6 an example of a passing particle with $\lambda =0.9$ and a trapped particle with $\lambda =0.99$, correspondingly. To observe the banana orbits characteristic of trapped particles in tokamaks, we transformed the VMEC equilibria to the BOOZ$\_$XFORM (Sanchez et al. 2000) format and calculated the particle trajectory for this field, where the angular coordinates are Boozer coordinates similar to pyQSC. The plots in Boozer coordinates in panels (d–f) of each figure show the comparison between the orbits in these two fields, which exhibit circular curves for passing orbits and bean or banana-shaped curves for trapped orbits. We also show a view from above of the orbits in the pyQSC field (qsc) and in the VMEC equilibrium (vmec). In this case, the same inputs for the particles result in the same type of orbits, but this does not have to be the case as in different aspect ratios, we have different values of initial $B$ and, therefore, for some values of $\lambda$, some orbits are trapped and some are not for different aspect ratios.

Figure 5. Comparison between alpha particle orbits in pyQSC (qsc), VMEC (vmec) and BOOZ$\_$XFORM (booz) QA equilibria with initial position $(s, \theta, \varphi ) = (0.25, 0.1, 0.1)$ in Boozer coordinates and $\lambda =0.8$. (a–c) Evolution of the radial position in time. (d–f) Normalized poloidal view of the orbit. (g–i) Top view of the orbit in Cartesian coordinates. (a,d,g) $A=6.8$. (b,e,h) $A=9.1$. (c, f,i) $A=13.6$.

Figure 6. Comparison between alpha particle orbits in pyQSC (qsc), VMEC (vmec) and BOOZ$\_$XFORM (booz) QA equilibria with initial position $(s, \theta, \varphi ) = (0.25, 0.1, 0.1)$ in Boozer coordinates and $\lambda =0.99$. (a–c) Evolution of the radial position in time. (d–f) Normalized poloidal view of the orbit. (g–i) Top view of the orbit in Cartesian coordinates. (a,d,g) $A=6.8$. (b,e,h) $A=9.1$. (c, f,i) $A=13.6$.

Figure 7. QH VMEC equilibria generated from the ‘2022 QH nfp4 well’ from pyQSC with different major radius scalings. Top panels are the 3-D versions of the equilibria and bottom panels are the contour plots of the magnetic fields on the angular VMEC coordinates. (a) Aspect ratio $A=6.8$. (b) $A=9.1$. (c) $A=13.6$.

For higher $\lambda$, orbits have some deviations from one another for all aspect ratios, but the most striking differences happen when we have trapped particles, as can be seen in figure 6, where the average radial positions deviate significantly for the lower aspect ratio which implies larger differences in all the points of view of the orbit. Although no loss of particles is shown in the figures, it is easy to imagine that, for wider banana orbits, we would lose a particle for all equilibria, at least for the higher aspect ratios, making the prompt loss fraction similar between these configurations, despite small differences in the magnetic field.

The fact that the orbits in the BOOZ$\_$XFORM (booz) and VMEC (vmec) fields are quite different from each other, figure 6 shows that slight changes on the equilibrium, which can arise from interpolation errors or resolution differences, can cause large differences in the radial excursion of a given particle with the same initial conditions when we are further away from QS or omnigenity. Although increasing the resolution could improve the convergence between fields, the required increase in toroidal and poloidal modes to achieve convergence would substantially affect the computational costs of the procedure. Furthermore, VMEC is known to have decreased accuracy near the magnetic axis, so some level of disagreement is expected there. However, even though the orbits are not aligned with each other, the impact of this effect in the confinement of these particles is not obvious as the misalignment in the form of time delays or small average deviations may not have a big impact on the fraction of lost particles. Such small deviations will be assessed in the next section.

4.2. Quasi-helical symmetry

The orbit widths obtained for the QH stellarator have, in general, smaller amplitudes, which matches the observation of Landreman & Paul (2022) that the radial excursions are expected to decrease with a higher value of $|\iota - N|$. The lower amplitudes, which can be observed in figures 8 and 9, appear to be accompanied by higher frequencies of oscillation. Additionally, passing orbits such as those in figure 8 exhibit proportionally higher deviations of the VMEC and BOOZ$\_$XFORM orbits from the pyQSC ones, although in absolute value, they tend to be smaller than their QA counterparts in figure 5.

Figure 8. Comparison between alpha particle orbits in pyQSC (qsc), VMEC (vmec) and BOOZ$\_$XFORM (booz) QH equilibria with initial position $(s, \theta, \varphi ) = (0.25, 0.1, 0.1)$ in Boozer coordinates and $\lambda =0.8$. (a–c) Evolution of the radial position in time. (d–f) Normalized poloidal view of the orbit. (g–i) Top view of the orbit in Cartesian coordinates. (a,d,g) $A=6.8$. (b,e,h) $A=9.1$. (c, f,i) $A=13.6$.

Figure 9. Comparison between alpha particle orbits in pyQSC (qsc), VMEC (vmec) and BOOZ$\_$XFORM (booz) QH equilibria with initial position $(s, \theta, \varphi ) = (0.25, 0.1, 0.1)$ in Boozer coordinates and $\lambda =0.99$. (a–c) Evolution of the radial position in time. (d–f) Normalized poloidal view of the orbit. (g–i) Top view of the orbit in Cartesian coordinates. (a,d,g) $A=6.8$. (b,e,h) $A=9.1$. (c, f,i) $A=13.6$.

For the case where $\lambda =0.99$, as depicted in figure 9, distinct behaviours emerge based on the aspect ratio. When considering a smaller aspect ratio, the orbit in the NAE magnetic field is a purely trapped particle, while the remaining orbits oscillate between trapped and passing states throughout the temporal evolution. In this scenario, the total radial amplitude of the transitioning trajectories exhibits only a marginal increase compared with the purely trapped orbits. For the medium aspect ratio, the particle in an NAE field is now a passing one, while the other particles continue to demonstrate the previously observed behaviour, with radial amplitudes exceeding twice that of the passing trajectory. For the case of a larger aspect ratio, as all particles are passing ones, only small deviations are observed between orbits. It is also interesting to note that the banana widths are smaller in the QH case than in the QA case, as expected.

Another important aspect is that no signs of trapped particles with an average outward radial drift, $\Delta \psi > 0$, were observed for the orbits in the QH NAE fields. This may in part be caused by the increased quasi-symmetry in these fields or an expression of the simplicity of the field, not allowing for locally trapped or detrapping particles. Passing–trapped transitions can be caused by the finite banana width, which is usually a small effect, and by the misalignment of field maxima, which can be attributed to a deviation from perfect quasi-symmetry and can cause stochastic diffusion losses of approximately $10\,\%$ for particles born at half-radius (Beidler et al. 2001). The latter effect is thus more relevant to the QH configurations, where the degree of quasi-symmetry achieved is smaller. We note that such deviation is expected to stem mainly from differences in $B$ as the same code, equations and algorithms were used. Despite this fact, the lower transport due to increased rotational transform, an important feature of most QH stellarators, as it provides better confinement of energetic particles for lower time frames, is partially captured by the NAE.

4.3. Trapped–passing separatrix

As the alignment of orbits for the different equilibria appears to be appreciably worse for trapped particles than for their passing counterparts, it becomes of interest to try to estimate the threshold between the two of them taking into account their initial coordinates. We attempt exactly that using the first-order NAE analytical formula for the magnetic field.

Defining the energy and magnetic moment as $E= m ( v_\perp ^2 + v_\parallel ^2) / 2$ and $\mu = {m v_\perp ^2}/{2 B(s, \theta, \phi )}$, respectively, and defining $\lambda = E_{\perp _i} / E = \mu B_i / E$, where the subscript $i$ stands for the initial value given for the simulation, we arrive at

(4.1)\begin{equation} v_\parallel ^2 = v^2 \left(1 - \lambda \frac{B(s, \theta, \phi)}{B_i}\right). \end{equation}

For a trapped particle, there is a given critical radial position $(s_c, \theta _c, \varphi _c)$ at which $v_\parallel = 0$, leading to $B(s_c, \theta _c, \varphi _c)) = B_c = { B_i}/{\lambda }$. With a first-order near-axis expansion, we can write $B_i = B_0 (1 + a_A \sqrt {s_i} \bar {\eta } \cos {\theta _i})$ and $B_c = B_0 (1 + a_A \sqrt {s_c} |\bar {\eta }| )$, which we take to be a poloidal maximum of $B$ for this estimation, which leads to $| \cos {\theta } | =1$. This assumption results in an overestimation of the value of $B_c$, but simplifies the analysis, leading to an estimation of $\lambda$ at which we have the passing–trapped separatrix, $\lambda _s$, depending only on the initial position of the particle and the radial position where the parallel velocity is null, $s_c$. Taking advantage of the fact that confined orbits perform an oscillatory motion, we may assume $s_c$ should be between $s_i - 2\Delta s$ and $s_i + 2\Delta s$, where $2 \Delta s$ is an estimation of the maximum radial peak-to-peak amplitude of orbits. This is done since the starting point can be in the peak or the valley of the oscillation. Finally, we arrive at the following interval for the transition value $\lambda _s$:

(4.2)\begin{equation} \frac{(1 + a_A \sqrt{s_i} \bar{\eta} \cos{\theta_i})}{(1 + a_A \sqrt{s_i + 2 \Delta s} |\bar{\eta}| )} \leq \lambda_s \leq \frac{(1 + a_A \sqrt{s_i} \bar{\eta} \cos{\theta_i})}{(1 + a_A \sqrt{s_i - 2 \Delta s} |\bar{\eta}| )}. \end{equation}

To demonstrate how effective this interval is at predicting if a particle is trapped or passing, we traced particles in the pyQSC QA and QH equilibria of aspect ratio $A=9.1$ with linearly spaced initial conditions $\theta$, $\varphi$ and $\lambda$ for three different radial positions $s_i=0.25$, $0.5$ and $0.75$. We then check if the particles have $v_\parallel =0$ at any point of their orbit (trapped) or not (passing), discarding all the lost particles classified as passing orbits. This last step was performed to avoid false negatives, as some orbits could be trapped particles that are lost before achieving $v_\parallel =0$. The interval in (4.2) is depicted in the form of three thresholds, two of them corresponding to the lower and higher bounds on $\lambda _s$ and the third one corresponding to $s_c=s_i$. For the QA stellarator, we take $2 \Delta s$ to be 0.4 and, for the QH, 0.1 as the maximum amplitudes in our examples appear to be of this order.

The obtained results can be seen in figures 10 and 11, where the threshold lines match the changes in the trapped area in a quite accurate way both in the QA and the QH stellarators for different radial initial positions. It is also interesting to note that the interval encompasses almost all overlapping dots even when the interval is of a smaller magnitude, such as in the QH case even though it is the product of a few assumptions. That said, the middle threshold can still be useful to have a single predictor of the type of orbit in question, especially in the QH case, where the amount of overlapping points is substantially smaller.

Figure 10. Comparison of estimated interval for trapped to passing separatrix with computed results for a QA stellarator. Trapped ( red) and passing ( blue) particles are plotted as a function of the initial Boozer poloidal coordinate $\theta$ and of $\lambda =\mu B_i/ E$ for three different initial radial positions $s_i$. High and low thresholds correspond to interval bounds, while mid threshold corresponds to an average separatrix : (a) $s_i=0.25$ ; (b) $s_i=0.5$ ;(c) $s_i=0.75$.

Figure 11. Comparison of estimated interval for trapped to passing separatrix with computed results for a QH stellarator. Trapped (red) and passing (blue) particles are plotted as a function of the initial Boozer poloidal coordinate $\theta$ and of $\lambda =\mu B_i/ E$ for three different initial radial positions $s_i$. High and low thresholds correspond to interval bounds, while mid threshold corresponds to an average separatrix: (a) $s_i=0.25$; (b) $s_i=0.5$; (c) $s_i=0.75$.

In future works, this interval could be used to select a priori only passing or trapped particles for analysis. This would be interesting to do because trapped particle orbits are generally unconfined and may exhibit chaotic behaviour, and could be more easily studied with this approach. Furthermore, these particles represent a substantial part of alpha-particle losses and so, computing loss fractions of trapped particles only would be more computationally efficient without considerable information losses. However, as the interval $\Delta s$ was found only a posteriori, it is still necessary to find a way of estimating this value. An attempt to do exactly this is found in the next subsection.

4.4. Radial oscillation amplitude

Not only is it of interest to have an estimate of the maximum oscillation amplitude $\Delta s$ to complete the above study, but also to have an expression for the dependence of the amplitude with $\lambda$ and $s_i$. To estimate these quantities, we take advantage of the fact that, for a general quasi-symmetric stellarator, the conserved momentum $p_\eta$ can be obtained from the Lagrangian in (3.3) through a coordinate transformation from $(s, \theta, \varphi )$ to $(s, \chi, \eta )$. This momentum can be written in the following way:

(4.3)\begin{equation} p_\eta = q \frac{N}{M} \psi - q \psi_p + \frac{m v_\parallel}{B(s, \theta, \varphi)} \left(G + \frac{N}{M}I\right), \end{equation}

and is the starting point of our estimation. Taking the first-order NAE, we approximate $I = 0$, $\psi _p=\int \iota \,{\rm d}\psi = \iota _0 \psi$ and $G = G_0= {L B_0}/(2 {\rm \pi})$, where $L=\int _0^{2{\rm \pi} }{\rm d}\phi \ell '$ is the axis length and $\ell$ is the arc length along the field line. Additionally assuming $\iota _N= \iota _0 - N$ and $M=1$, and using (4.1), we find an expression for the radial excursion of the particle in flux surfaces

(4.4)\begin{equation} \psi ={-} \frac{p_\eta}{q \iota_{N}} \pm \frac{m v L B_0}{2 {\rm \pi}q \iota_{N} B(s, \theta)}\sqrt{1 -\frac{\lambda B(s, \theta)}{B_i}}, \end{equation}

which depends on the magnetic field $B$, which in turn depends on the radial position $s$ and poloidal angle $\theta$, an implicit equation that would have to be computationally solved. However, we can estimate the amplitude of the orbits by applying the linearization $\Delta \psi / \Delta B = \partial _B \psi$ to (4.4) and taking $\partial _B p_\eta = 0$. For the estimation of the variation of $B$, we take into account that it depends on $s$ and $\theta$, and thus sum its variation in both coordinates

(4.5)\begin{align} \Delta B & = \int_0^{\rm \pi} \partial_\theta B \,{\rm d}\theta + \int_{s_i}^{s_{i}+ \Delta s_{\max}} \partial_s B \,{\rm d} s \nonumber\\ & = (B(s_i, {\rm \pi}) - B(s_i, 0)) + (B(s_i + \Delta {s_{{\rm avg}}}, \theta_i) -B(s_i, \theta_i)) \nonumber\\ & = |\bar{\eta}| a_A B_0 (2 \sqrt{s_i} + (\sqrt{s_i + \Delta {s_{\max}}} - \sqrt{s_i}) \cos{\theta_i}). \end{align}

Using the expression for $\Delta B$ in (4.5), applying the partial derivative of $B$ and making $B(s,\theta ) = B_0$, we obtain an estimation of the amplitude of oscillation of the flux surface coordinate $\psi$ which, written as the normalized radial coordinate $s=\psi / \psi _b$, with $\psi _b = B_0 a_A^2 / 2$, becomes

(4.6)\begin{equation} \Delta s = \frac{ m v L \bar{\eta}}{{\rm \pi} q \iota_{N_0} a_A B_0} \frac{1 - \lambda B_0/ (2 B_i)}{\sqrt{1-\lambda B_0/B_i}} ( 2 \sqrt{s_i} + (\sqrt{s_i + \Delta {s_{{\rm avg}}}} - \sqrt{s_i}) \cos{\theta_i}). \end{equation}

This expression leads to results that approximate the computed ones up until the passing–trapped separatrix, as can be seen in figures 12 and 13, starting to deviate from it afterwards. This appears to be the limit of the applied linearization. To extend the region of applicability of our estimation, it may be possible to go to higher orders on the approximation. The expression appears to capture the differences in amplitude for different stellarator configurations and different initial conditions.

Figure 12. Comparison of computed with estimated orbital radial amplitude $\Delta s$ in a QA configuration as a function of $\lambda = \mu B_i / E$ for different initial poloidal and radial coordinates. Vertical lines indicate the average passing–trapped separatrix for each $\theta _i$, where the two quantities start to diverge : (a) $s_i=0.25$ ; (b) $s_i=0.5$ ; (c) $s_i=0.75$.

Figure 13. Comparison of computed with estimated orbital radial amplitude $\Delta s$ in a QH configuration as a function of $\lambda = \mu B_i / E$ for different initial poloidal and radial coordinates. Vertical lines indicate the average passing–trapped separatrix for each $\theta _i$, where the two quantities start to diverge : (a) $s_i=0.25$ ; (b) $s_i=0.5$ ; (c) $s_i=0.75$.

Concerning a reasonable estimation of the maximum of the radial amplitude of the motion, it would be necessary to neglect the $\theta$ component of $\Delta B$. However, the maximum value would remain outside this regime so it may not be prudent to use the above expression to obtain the wanted value directly. The value of the proposed expression at the separatrix could instead give an initial guess of the appropriate value.

In conclusion, through examining single-particle orbits in different equilibria, we have gained critical insights into the dynamics of particles in NAE and ideal MHD magnetic fields. As we move forwards to the analysis of particle ensemble behaviour and loss fractions, we expect to improve our understanding of the collective behaviour of particles, providing a more holistic view of the confinement process.

5. Particle ensembles

We now change the focus from single-particle trajectories to high-energy particle ensembles. The analysis of these groups gives us a more comprehensive look at how well particles are confined and can reveal patterns that might not be evident when examining single-particle dynamics. In this part of our study, we focus on calculating loss fractions, which are the percentage of particles that leave the magnetic equilibrium over a certain time frame. These are crucial for assessing how effective the magnetic fields are at keeping particles in place, making them a helpful metric to inform improvements in our confinement techniques.

To determine loss fractions, large groups of particles are initialized for a given duration and the percentage of these particles that cross the last closed flux surface is accounted for. It is important to acknowledge that the real values of loss fractions in an actual device are multifactorial, and the computations performed in this work only take into account some of the relevant physical factors. For a thorough analysis of non-collisional loss fractions, simulations should account for a wide range of random starting conditions and incorporate a density distribution for alpha particles. A realistic account of losses in a given configuration should also take into account that not all of the particles crossing the set boundary would be lost, as some could later reenter the plasma, in addition to taking into account the effects of collisions and turbulence. The former effect was estimated to affect the losses up to ${\sim }10\,\%$ in the LHD (Miyazawa et al. 2014). However, for a comparison between the losses in different fields, we are less concerned with the realistic global loss values than we are with having similar results between the obtained values, so we compute the loss fractions for different initial radial positions, and we have linearly spaced initial poloidal and toroidal VMEC coordinates and $\lambda$ in the following ranges: $[0, 2 {\rm \pi}]$, $[0, 2 {\rm \pi}/n_{fp}]$ and $[0, 1]$, respectively. For every set of these initial conditions, $s_{v_\parallel }$ is set as $-$1 and $+$1, so our number of particles $n_p$ is equal to two times the product of $n_\theta$, $n_\phi$ and $n_\lambda$, where $n_x$ is the number of different possible initial values of the initial quantity $x$. To ensure the initial positions matched in the VMEC and pyQSC equilibria, a root finder was designed to compute the corresponding $\phi _0$ for every $\phi$, with $\varphi$ calculated subsequently. The remaining NAE coordinates were calculated following the reverse of the procedure described in § 2. This strategy was employed to minimize the arguments passed to the VMEC loss fraction algorithm, as this was the computational bottleneck in the workflow. Additionally, as discussed above, the time frame where the NAE seems to be more relevant is the prompt loss one, and so we will simulate the orbits up until $5 \times 10^{-3}$ s. We note that for all the results obtained in this work, the computation of loss fractions is two orders of magnitude faster for the pyQSC magnetic fields when compared with the VMEC calculations, mainly due to the simple closed form solutions available of the first fields.

In the forthcoming discussions, we will connect our earlier findings on individual particle behaviour to this larger-scale evaluation. By doing so, we aim to better understand the level of particle confinement in MHD and NAE systems, and their differences.

5.1. Quasi-axisymmetry

For the calculation of the loss fractions on a QA stellarator, we recur to the same baseline equilibrium precise QA, as in the last section, with an aspect ratio of $A= 9.1$. We then trace ensembles of particles with initial radial positions $s_i=0.25$, $0.5$ and $0.75$, facilitating an investigation into the dependence of result correspondence upon the radial distance from the equilibrium axis.

In figure 14, we present the results of the aforementioned ensembles, where we can see an overall agreement over all initial radial positions. The fact that this agreement is present despite differences in the radial oscillations of particles may be attributed to the averaging out of VMEC orbits that go slightly above and below the NAE orbits for different values of $\phi$. This conclusion is supported by the fact that increasing the value of $n_\phi$ improves the agreement of the loss fraction.

Figure 14. Comparison between alpha particle ensemble loss fractions in pyQSC (NA) and VMEC (VMEC) with $n_p=9000$, and linearly spaced initial conditions $\theta _i$ and $\phi _i$ in VMEC coordinates as well as $\lambda$ for the aspect ratio $A=9.1$ QA equilibria with initial radial positions (a) $s_i=0.25$, (b) $s_i=0.5$ and (c) $s_i=0.75$.

There is an uncanny agreement for all presented flux surfaces up until the time mark of $t=5 \times 10^{-4}$ s, where the values for the two equilibria start to diverge. This is most likely where other effects start being relevant for this stellarator other than first-orbit losses. The relative difference in loss fractions is the highest for $s_i=0.25$, where the pyQSC equilibrium has no prompt losses, but the percentage on the VMEC losses stays below $1\,\%$ on the analysed time frame. For the ensemble most distant from the axis in figure 14, the amount of lost particles appears to match for the entire measured time, indicating the loss cone particles are the most important loss mechanism at that flux surface.

These results seem to indicate that the computation of loss fractions for initial time intervals could be performed with NAE fields with a low loss in performance and a reduction of the computation time of two orders of magnitude. Each of the $n_p=9000$ ensembles took approximately four hours to compute in the VMEC field and two minutes in an openMP-only parallelized code running on eight cores.

5.2. Quasi-helical symmetry

A similar procedure as in the preceding subsection was followed for the 9.1 aspect ratio QH stellarator used in the single-particle tracing. As the orbit amplitude is significantly narrower in this equilibrium, first-orbit losses are only seen in its outer regions. For this reason, we present only the results of the ensembles with initial radial positions up from $s_i=0.5$, as results for inner radial starting points do not vary qualitatively, only quantitatively. As the frequency of oscillation in this magnetic field was higher, the divergence between equilibria arose sooner too. Therefore, the integration only goes up to $t=5 \times 10^{-4}$ s.

Although, as shown in figure 15, the behaviour of the first few orbits is well captured by the NAE, in the present case, there are no first-orbit losses and, therefore, the loss fraction remains zero until $t=1 \times 10^{-4}$ s for the cases of $s_i=0.5$ and $0.75$, and later diverges when particles start exhibiting an average radial drift outwards in the VMEC equilibrium. In the case of the ensemble generated with initial radial position $s_i=0.9$, there is a region of time where the losses in the VMEC magnetic field are emulated by the ones in the pyQSC field. However, they start to diverge before the $t=1 \times 10^{-4}$ mark. The fact that the convergence is limited in time and space makes the argument for optimization in this case more limited. However, the fact that the QH VMEC equilibrium exhibited a larger deviation from quasi-symmetry may also be responsible for the accelerated divergence between the loss fraction results, as non-vanishing average radial drifts can be a consequence of local loss of symmetry. It can also be a result of particles transitioning from the passing to the trapped regime in VMEC without doing so in pyQSC fields such as in figure 9.

Figure 15. Comparison between alpha particle ensemble loss fractions in pyQSC (NA) and VMEC (VMEC) with $n_p=9000$, and linearly spaced initial conditions $\theta _i$ and $\phi _i$ in VMEC coordinates as well as $\lambda$ for the aspect ratio $A=9.1$ QH equilibria with initial radial positions (a) $s_i=0.5$, (b) $s_i=0.75$ and (c) $s_i=0.9$.

Although the effectiveness of the tool explored throughout this work appears to be limited to prompt losses, it is important to relate our findings to the results of LeViness et al. (2022) and Velasco et al. (2023). The fact that these results point to some correlation between prompt losses and proxies for losses suggests that addressing the mechanisms responsible for these kinds of losses may also lead to a reduction of other loss mechanisms that happen at later times. This way, increasing the efficiency for prompt losses optimization can be beneficial beyond the initially expected time frames. Additionally, we note that the application of the followed methodology is more insightful than simply pursuing QS in a stellarator, providing a direct loss measurement, and more accurate than the aforementioned proxies for initial losses as it includes finite banana width effects. Although lacking in accuracy when compared with direct GC simulations in MHD equilibria (Paul et al. 2022), it benefits from increased computational speed as expressed in this work.

6. Conclusions

In this work, we started by comparing GC orbits computed by the Euler–Lagrange equations of motion implemented in the gyromotion library gyronimo with those obtained by the symplectic formulation implemented in the loss-fraction assessment code SIMPLE. The explicit methods in gyronimo compared well with the implicit symplectic one, showing good levels of numerical agreement. Subsequently, single-particle tracing was performed in different aspect ratio reactor scale stellarators. The discrepancy between the NAE fields and the MHD ones such as the VMEC and BOOZ$\_$XFORM did not have a noticeable first-order effect on the passing particles independently of the aspect ratio and stellarator, with growing deviations for lower aspect ratios especially in the QH equilibrium. More noticeable disparities were found for the trapped particles, where passing–trapped transitioning particles, radially drifting particles and other effects were observed for the MHD equilibria in contrast with the NAE ones, where no such phenomena were identified. Once more, these effects were more prominent in the QH magnetic fields. These discrepancies were expected due to the decreased levels of quasi-symmetry in the MHD fields. Additionally, effective expressions for estimating the radial amplitude of particles’ orbits and the passing–trapped separatrix were found for both QA and QH stellarators in the NAE.

Following the tendencies of the single-particle tracing orbits, for each equilibrium studied, the loss fractions of large groups of particles at different initial radial positions agree only for losses of orbits whose radial amplitude is larger than their distance to the edge of the magnetic field. For the QA stellarator, this agreement lasts up until $t=5 \times 10^{-4}$, which could be used to do a first optimization of stellarators or at least eliminate classes of stellarators with too many first-orbit losses. Unfortunately, for the QH stellarator, these kinds of losses do not seem to be relevant for the majority of initial positions, as radial amplitudes tend to be smaller, which is concordant with the expectation for a larger value of $\iota$. It is important to emphasize the loss fraction results are obtained up to two orders of magnitude faster for the NAE particle tracing. Although faster performance could be attained in the different tracings by specializing the code to a given field, this is outside the scope of this work.

Future work with loss fractions using the NAE should encompass a larger diversity of equilibria, including equilibria generated from relevant MHD fields by fitting them to near-axis fields to compare with established results. The implementation of the loss fraction estimations in the optimization of directly constructed equilibria could also be studied as a way of producing a good initial configuration for conventional optimization.