2. The near-axis expansion

In this section, we follow the notation of Landreman & Sengupta (Reference Landreman and Sengupta2019) and introduce the near-axis expansion using the inverse approach first presented by Garren & Boozer (Reference Garren and Boozer1991a). 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 Reference Landreman and Sengupta2018). In quasi-symmetry, $\boldsymbol {r}$ can be written as (Landreman & Sengupta Reference Landreman and Sengupta2019)

(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 (Reference Landreman and Sengupta2019) 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 Reference Landreman2021) and with VMEC (Hirshman Reference Hirshman1983), 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 (Reference Landreman and Sengupta2018). 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. Reference Bader, Drevlak, Anderson, Faber, Hegna, Likin, Schmitt and Talmadge2019a; Henneberg et al. Reference Henneberg, Drevlak, Nührenberg, Beidler, Turkin, Loizu and Helander2019), 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 Reference Albert, Kasilov and Kernbichler2020a). 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 Reference Cary and Brizard2009):

(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 Reference Littlejohn1983):

(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 Reference Cary and Brizard2009). 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 Reference Landreman and Paul2022).

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. Reference Albert, Kasilov and Kernbichler2020b), 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 (Reference Albert, Kasilov and Kernbichler2020c), 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. (Reference Albert, Kasilov and Kernbichler2020a).

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 Reference Ahnert and Mulansky2011), 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. Reference Najmabadi and Raffray2008), 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 Reference Jorge and Landreman2021). 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 (Reference Landreman and Paul2022) 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. Reference Albert, Kasilov and Kernbichler2020c). 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 (Reference Landreman and Paul2022) and the one for QH stellarators was the four field period QH vacuum configuration with magnetic well from § 5.4 of Landreman (Reference Landreman2022). 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. (Reference Landreman, Sengupta and Plunk2019). 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. Reference Sanchez, Hirshman, Ware, Berry and Spong2000) 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 (Reference Landreman and Paul2022) 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. Reference Beidler, Kolesnichenko, Marchenko, Sidorenko and Wobig2001). 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. Reference Miyazawa, Suzuki, Satake, Seki, Masaoka, Murakami, Yokoyama, Narushima, Nunami and Goto2014). 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. (Reference LeViness, Schmitt, Lazerson, Bader, Faber, Hammond and Gates2022) and Velasco et al. (Reference Velasco, Calvo, Sánchez and Parra2023). 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. Reference Paul, Bhattacharjee, Landreman, Alex, Velasco and Nies2022), it benefits from increased computational speed as expressed in this work.