1. Introduction
Minimizing the losses of fusion alpha particles is a major goal in the optimization of stellarator magnetic confinement configurations for self-heated fusion plasmas. Recent results (Goodman et al. Reference Goodman, Mata, Henneberg, Jorge, Landreman, Plunk, Smith, Mackenbach and Helander2022; Landreman, Buller & Drevlak Reference Landreman, Buller and Drevlak2022; Landreman & Paul Reference Landreman and Paul2022; Paul et al. Reference Paul, Bhattacharjee, Landreman, Alex, Velasco and Nies2022) demonstrate that it is possible to achieve excellent fusion alpha confinement in stellarator reactor configurations optimized towards precise quasi-symmetry or quasi-isodynamicity. Retaining this quality of alpha confinement while also fulfilling other objectives such as magnetohydrodynamic stability, low turbulent transport or coil complexity is a major challenge. Therefore, criteria for good alpha confinement that are more general than quasi-symmetry and quasi-isodynamicity are still of high interest. Here, we present and assess such criteria.
One way to obtain alpha loss fractions in a stellarator configuration consists of directly tracing guiding-centre orbits and (optionally) taking into account collisions via Monte Carlo operators. The slowing-down time of fusion alphas is in the range of 0.1–0.5 seconds. To decide whether a single guiding-centre orbit remains confined during this time typically requires tracing for ${>}10^4$
bounce periods. Numerical estimation of loss fractions from a statistical ensemble of $10^4$
particles and hundreds of time steps per period therefore requires billions of magnetic field evaluations. Despite recent acceleration by symplectic integration (Albert, Kasilov & Kernbichler Reference Albert, Kasilov and Kernbichler2020b) and Poincaré plot classification (Albert, Kasilov & Kernbichler Reference Albert, Kasilov and Kernbichler2020a; Kamath Reference Kamath2022), direct computation of alpha losses remains computationally expensive compared with other components of a stellarator optimization loop. While there has been recent progress on accelerating prompt loss computations (Velasco et al. Reference Velasco, Calvo, Mulas, Sánchez, Parra and Cappa2021), there remains room for improvement of optimizing also late and collisional alpha energy losses. Since pitch-angle scattering becomes relevant only after alpha particles have slowed down significantly (see, e.g. Ho & Kulsrud Reference Ho and Kulsrud1986), collisionless proxies can be useful also for this purpose.
To provide more efficient and accurate proxies, two novel classifiers are introduced and their performance and reliability first assessed on a test case in the optimized stellarator configuration of Drevlak et al. (Reference Drevlak, Beidler, Geiger, Helander and Turkin2014). In the terminology of statistics, collisionless orbit classification checks the hypothesis of whether a guiding-centre orbit is lost by first testing the weaker hypothesis whether it is chaotic (potentially lost). This is done in a way analogous to a screening test for a disease. If a classifier wrongly identifies a regular orbit (never lost) as chaotic, this is known as a false positive or type I error. Another more accurate (and expensive) test will subsequently be required to come to a final conclusion. In our case, we directly check whether one of the pre-screened orbits deemed chaotic will be finally lost. More type I errors produce more (computational) cost, but the final result remains unchanged. In contrast, a wrong identification of a chaotic orbit as regular, is known as a false negative or type II error. In that case, the orbit is prematurely counted as confined, even though it may be lost later. While this kind of error must be avoided for exact final results, it may also be permissible as long as the error rate is sufficiently small, and the final error does not exceed statistical fluctuations. We compare two classifiers with the existing Minkowski classifier in terms of these error modes. The first classifier is based on topological properties of a subset of regular orbits in phase space, and the other one on the variation of the parallel adiabatic invariant.
New classifiers are presented in § 2. Based on the increased computational performance of these classifiers, it becomes possible to analyse the regular and chaotic regions with a fine resolution across the whole phase space of collisionless fusion alpha orbits in a stellarator configuration in comparably short time. Such analyses are presented in § 3 for five optimized stellarator reactor configurations (Subbotin et al. Reference Subbotin, Mikhailov, Shafranov, Isaev, Nührenberg, Nührenberg, Zille, Nemov, Kasilov and Kalyuzhnyj2006; Drevlak et al. Reference Drevlak, Beidler, Geiger, Helander and Turkin2014, Reference Drevlak, Beidler, Geiger, Helander, Henneberg, Nührenberg and Turkin2018; Henneberg et al. Reference Henneberg, Drevlak, Nührenberg, Beidler, Turkin, Loizu and Helander2019; Landreman & Paul Reference Landreman and Paul2022). The results are compared with direct computations of collisionless and collisional alpha particle and energy loss fractions from the guiding-centre code SIMPLE (Albert et al. Reference Albert, Kasilov and Kernbichler2020a,Reference Albert, Kasilov and Kernbichlerb). Based on the quasi-axisymmetric case, the role of orbital resonances is pointed out, with more details described in a recent paper (Albert et al. Reference Albert, Rath, Buchholz, Kasilov and Kernbichler2022).
2. Methods
Here, we build on the idea of classification of collisionless orbits with help of Poincaré plots. In an earlier work (Albert et al. Reference Albert, Kasilov and Kernbichler2020a), distinguishing regular and chaotic orbits has been based on box counting as an estimate of the Minkowski fractal dimension. This method required approximately 10 000 footprints in a Poincaré plot to work reliably. Our aim is to devise alternative algorithms which will distinguish regular trapped orbits from chaotic ones using a much lower number of footprints.
We start from ideal orbits, which approximately follow poloidally closed contours of the parallel adiabatic invariant $J_\parallel$
in quasi-isodynamic (QI) configurations (Gori, Lotz & Nührenberg Reference Gori, Lotz and Nührenberg1996). We call them ‘ideal’ because all orbits should be like this in a perfect QI configuration. The same definition of ideal orbits applies also to the quasi-helical (QH) configuration but not to the quasi-axisymmetric (QA) configuration where the symmetry direction is toroidal and, therefore, contours must be closed toroidally, not poloidally. Up to this change, the results for QI configurations apply also to QA configurations. Regular orbits are a superset of ideal orbits, as their footprints are required to lie on closed curves which are not necessarily single-valued functions of the relevant angle.
An example Poincaré plot of the banana tips (The Poincaré cut is a hypersurface of parallel velocity $v_\parallel =0$
with $v_\parallel$
changing sign from negative to positive) of an ideal orbit in QI configuration is shown in figure 1 in the $(\theta,s)$
plane, where $s$
is a flux surface label and $\theta$
is a poloidal angle. If $\theta$
is an angle in quasi-symmetry direction (toroidal angle in case of QA configurations), and quasi-symmetry is perfect, banana tips of any orbit lie on a curve $s(\theta )=\mathrm {const}$
and therefore all orbits are ideal. Also in this figure, the lowest-order parallel adiabatic invariant $J_\parallel$
is shown that has been computed in two ways as follows:
Figure 1. Poincaré plot of ${\approx }10.000$
footprints in $(\theta,s)$
plane (a) and plot of $(\theta,J_\parallel )$
pairs for an ‘ideal’ orbit of a 3.5 MeV alpha particle in the QI stellarator reactor configuration of Drevlak et al. (Reference Drevlak, Beidler, Geiger, Helander and Turkin2014). First return invariant, $J_\parallel ^r$
, and bounce invariant $J_\parallel ^b$
are shown with black and grey dots, respectively.
The ‘first return’ invariant $J_\parallel ^r$
is computed along the real orbit with finite Larmor radius (FLR) as an integral of the squared parallel velocity $v_\parallel ^2$
over the return time to the Poincaré cut $\tau _r$
i.e. along the orbit connecting two subsequent footprints. The ‘bounce’ invariant $J_\parallel ^b$
is computed the usual way, along the field line, i.e. in the limit of zero Larmor radius $\rho _L$
. In the latter approximation, the orbit is closed within a bounce time $\tau _b$
which differs from $\tau _r$
by first-order terms in Larmor radius, unless the FLR orbit undergoes a transition. It can be seen that our ‘ideal’ orbit actually does not follow exactly the contours of either of these two. Although the changes in lowest-order $J_\parallel$
are within the fourth digit, they are of the same order for both definitions. Actually, this variation is because lowest-order $J_\parallel$
differs from the exact invariant by the first-order terms in Larmor radius (Hastie, Taylor & Haas Reference Hastie, Taylor and Haas1967). We will use this property later to introduce a classification method based on the relative variation in $J_\parallel$
.
Despite small variations in $J_\parallel$
, it is clearly seen in figure 1 that all the footprints are located on a smooth curve $s=s(\theta )$
– the cross-section of the Kolmogorov-Arnold-Moser (KAM) surface (‘drift surface’) by the Poincaré cut $v_\parallel =0$
. This curve is an invariant manifold of the Poincaré map, each mapping maps this curve onto itself, but the points get displaced along the curve. Let us use this property for orbit classification. According to Poincaré recurrence theorem, an orbit performing motion within a closed volume returns to any small vicinity of the starting point after a finite time, $\tau _R$
, which depends on the size of this vicinity, and it does this subsequently an infinite number of times. Besides the vicinity size, the recurrence time $\tau _R$
depends also on the type of the orbit. For the ideal regular orbits such as the one shown above it has a near-periodic property which we discuss and use below.
Figure 2 shows a comparison of different types of non-ideal orbits. The orbits on the left side (nonlinearly trapped in a resonance and transient regular orbit) are regular and therefore absolutely confined despite being non-ideal. The orbit trapped deeply in a drift-orbit resonance (see § 3.1) is similar to the non-ideal orbits shown in figures 4 and 5. Such orbits appear as narrow island chains and are thus prone to mis-classification due to their similarity to ideal orbits and relatively low variations in the parallel invariant. In contrast, the footprints of regular transient orbits, i.e. orbits which undergo trapping class transitions of the second kind during their precession (see the next paragraph for more details), switch between distinct islands in phase space. Since trapped resonant orbits and transient regular orbits are both regular, mis-classification does not affect final results for losses. Other non-ideal orbits (near separatrix of resonance and transient chaotic) are potentially lost and have to be traced until their slowing down time in calculations of loss statistics.
Figure 2. Examples of ‘non-ideal’ orbits for 3.5 MeV alphas in the QI configuration of Drevlak et al. (Reference Drevlak, Beidler, Geiger, Helander and Turkin2014) started from flux surface $s=0.6$
. Orbits nonlinearly trapped in a resonance (a) and some transient orbits (c) are regular and therefore confined. In contrast, orbits close to a resonance separatrix (b) and transient chaotic orbits (d) are possibly lost.
Here, we distinguish trapping classes by the number of local field minima traversed over the first return time $\tau _r$
. With this definition, there are two kinds of trapping class transitions possible during the precession. Within a transition of the first kind, minima are split (or merged) at the bottom of the magnetic well so that the normalized perpendicular kinetic energy, ${\mathcal {E}}_\perp =v_\perp ^2/v^2$
, is smaller than the one at newly appearing (or disappearing) local field maxima in the domain of bounce (first return) motion. These transitions have no effect on the orbits and change the trapping class only formally. Within a forward transition of the second kind (‘retrapping’), ${\mathcal {E}}_\perp$
reaches unity at some local maximum in the domain of bounce motion thus splitting it into two domains with only one chosen then by the particle in accordance with its bounce phase. This bounce phase $\theta ^2=2{\rm \pi} \tau /\tau _r$
parametrizes the particle position along the magnetic field line, where $\tau$
is the time elapsed after the previous return (see, e.g. Albert et al. Reference Albert, Heyn, Kapper, Kasilov, Kernbichler and Martitsch2016). The backward transition (‘detrapping’) is the inverse process where the domain of bounce motion abruptly increases. In both transitions of the second kind, the guiding-centre orbit crosses the separatrix of one-dimensional bounce motion (Kolesnichenko, Lutsenko & Tykhyy Reference Kolesnichenko, Lutsenko and Tykhyy2022) which results in the abrupt change of the parallel adiabatic invariant. Thus, ‘transient’ orbits are defined as those which undergo these class transitions of the second kind.
2.1. Topological orbit classifier
We now want to develop the first idea used to classify orbits based on the already mentioned near-periodic property of the recurrence time $\tau _R$
for ideal orbits. Let us first enumerate subsequent footprints $(\theta _k,s_k)$
starting from $k=0$
and fix the interval $\theta _0<\theta <\theta _1$
between the first two footprints. An orbit will return to this interval after $N_1$
mappings, where $N_1$
is the first recurrence number. Subsequent recurrence numbers $N_j$
are defined such that $\theta _0<\theta _{N_j}<\theta _1$
. For an ‘ideal’ orbit, all $N_j$
are inside the interval $N_1-1\le N_j-N_{j-1} \le N_1$
. The proof is apparent from figure 3. The segment of the invariant line between the initial two footprints 0 and 1 is mapped onto the segment between footprints $N_1-1$
and $N_1$
after $N_1-1$
Poincaré mappings. By a continuity argument, there exists a point $X$
within the original segment $[0:1]$
which is mapped exactly onto point 0 after $N_1-1$
mappings. Therefore, all points of the segment $[0:1]$
located between points $X$
and 1 are mapped onto the sub-segment between points $0$
and $N_1$
performing a recurrence after $N_1-1$
mappings. Points located between points 0 and $X$
are mapped onto the sub-segment between points $N_1-1$
and 0. The latter points require one more mapping in order to arrive in a sub-segment between points $N_1$
and 1. These points perform a recurrence after $N_1$
mappings. Thus, if $N_{j-1}$
is a recurrence number $j-1$
such that $\theta _0<\theta _{N_{j-1}}<\theta _1$
, the next recurrence number will be either $N_j=N_{j-1}+N_1-1$
or $N_j=N_{j-1}+N_1$
, so that both cases satisfy $N_1-1\le N_j-N_{j-1} \le N_1$
with $j\ge 2$
.
Figure 3. First two (red) and last two (blue) footprints of the mapping completing the first recurrence. The point $X$
is plotted in magenta.
By a similar argument, each second recurrence satisfies $N_2-1\le N_{2k}-N_{2(k-1)} \le N_2$
with $k\ge 2$
, and, more generally, each $n$
-th recurrence satisfies
with $k\ge 2$
. Note that a recurrence $N_j$
approximately corresponds to $j$
poloidal precession turns (the higher $j$
, the more accurate is this relation). Therefore, prompt orbit losses occur within the first recurrence.
The simplest algorithm to detect ideal orbits using recurrence (2.2) for $n=1$
is then straightforward. Every orbit is followed over $M$
precession turns, and after each turn $j$
a recurrence number $N_j$
is checked to be within $N_1-1\le N_j-N_{j-1} \le N_1$
. If this condition is fulfilled for all $M$
precession turns, an orbit is classified as ‘ideal’. Otherwise, it is classified as ‘non-ideal’. The results of this classification are compared with classification using Minkowski dimension in table 1. The table compares the classification of the orbits for different numbers of precession turns.
Table 1. Comparison of classification by Minkowski dimension and by topological classification using a simple recurrence relation.
As not all regular orbits are ideal, but all ideal orbits are regular, the topological classifier may correctly identify a regular orbit as non-ideal (figure 2). This does not lead to errors in the directly computed collisionless alpha loss result using the topological classifier for speed up, but only increases computation time, as non-ideal/chaotic orbits are traced to the end. In contrast, an orbit wrongly identified as ideal by the topological classifier leads to an underestimation of losses.
In the absence of errors, all chaotic orbits should be identified as non-ideal, however, this is not immediately the case. It can be seen that chaotic orbits wrongly counted as ideal by the topological classifier disappear after 16–32 precession turns. However, depending on the number of precession turns, some non-ideal, regular orbits are identified as ideal. While this does not affect prediction of losses, it should be kept in mind for applications where this distinction is important.
The single remaining mismatch after 32 turns in table 1 has been manually identified as a simultaneous mis-classification of a regular non-ideal orbit due to a type I error by the Minkowski method and a type II error by the topological classifier, see figure 4. With larger sample sizes, an accurate classifier might require 128 turns. We have used this number here to demonstrate how the classification results stabilize when going towards a ‘large’ number of footprints. The number of footprints for one turn is approximately 40, therefore one may require up to 4000 footprints – not much less than for the Minkowski classifier. If one decides to tolerate some mis-classification, note that for 8 turns, there are still 3 unidentified orbits (error around 1 %). As expected, some orbits classified as regular by Minkowski dimension are ‘non-ideal’ in the recurrence classification. Those orbits cannot be ‘ideal’ orbits because ideal orbits in this algorithm can be misclassified only due to numerical errors which are at the level of computer accuracy in our case.
Figure 4. Poincaré plots in $(\theta,s)$
(a,c) and $(\theta,J_\parallel )$
(b,d) planes for two regular orbits classified as ‘ideal’ by condition (2.3). A typical ideal and therefore regular orbit is shown in the upper panel, and a type I mis-classification by Minkowski dimension of a regular but non-ideal orbit is shown in the lower panel.
The simplest topological classification algorithm may be improved in two ways. First, by checking relations (2.2) for all possible $n$
values which are allowed by the condition $2n< M$
, so that recurrence with the maximum $n$
value is checked at least once. Second, one may additionally check the monotonicity of the ordered footprint sequence described below.
Let us order the recurrence numbers $N_1, N_2, \ldots, N_M$
in increasing sequence with respect to $\theta$
. Namely, we order their indices $j=j_k$
where $k=1,2,\ldots,M$
so that $\theta _0<\theta _{N_{j_1}}<\theta _{N_{j_2}}<\cdots <\theta _{N_{j_M}}<\theta _1$
. For an ideal orbit, this sequence is kept in all intervals between the subsequent footprints,
where $k \ge 1$
, and we assume no periodic boundaries within intervals $[\theta _k, \theta _{k+1}]$
. Condition (2.3) obviously follows from mapping of the continuous segment $[\theta _0, \theta _{1}]$
onto the segment $[\theta _k, \theta _{k+1}]$
after $k$
mappings in case the cross-section of the KAM surface is a simply connected line defined by some single-valued function $s=s(\theta )$
. We apply this additional restriction to the orbits which have been qualified as ‘ideal’ by the improved (multiple) recurrence classification, in order to remove further type II errors.
The result of additional classification using the monotonicity condition (2.3) is compared with classification by Minkowski dimension in table 2. It can be seen that all chaotic orbits are sorted out by this criterion already after 8 precession turns (the remaining orbit is, as we remember, a mis-classification by the Minkowski dimension). A manual study of the first 47 orbits classified by monotonicity condition as ‘ideal’ shows that indeed, most of them are ideal, like the first orbit in figure 4. There is only one exception of the type shown in the lower panel of figure 4. This orbit is re-classified as non-ideal after 32 precession turns.
Table 2. Comparison of classification by Minkowski dimension and by using the topological classifier with multiple recurrence and monotonicity criteria.
The error mode is shown in figure 5 where an orbit close to the periodic orbit (invariant axis) is plotted for 16, 32 and a much larger number of precession turns resulting in 609, 1217 and 8918 footprints, respectively. With 16 precession turns, the orbit is wrongly identified as ideal by both, recurrence and monotonicity conditions. For 32 turns, the version with multiple recurrence and monotonicity (2.3), becomes correct, while using simple recurrence alone still leads to a mis-classification.
Figure 5. Example of a regular non-ideal orbit that is mis-classified as ideal by the topological classifier (zoomed view to the right). First 609, 1217 and 8918 footprints are plotted with red circles, blue and black dots, respectively.
2.2. Parallel invariant classifier
A different classification method relies on the approximate conservation of the parallel adiabatic invariant $J_\parallel =J_\parallel ^r$
. As we have seen, for the ideal orbits and island-type orbits resulting from nonlinear trapping into the drift-orbit resonance, variation of $J_\parallel$
is rather small, below $1\,\%$
. More generally, also for most of non-ideal but regular orbits, the relative variations in $J_\parallel$
remain small. Exceptions are regular transient orbits (figure 2c) which are rather rare. This is why this classifier may be used to distinguish regular from chaotic orbits, rather than ideal from non-ideal ones.
Here, we have implemented classification with lowest-order $J_\parallel$
via a condition of regular orbits
where $J_\parallel ^{(k)}$
are the values of $J_\parallel ^r$
defined by (2.1a) for the first return period starting from footprint number $k$
. The result of classification using the condition (2.4) with ${\rm \Delta} J_\parallel$
chosen as $1\,\%$
of its initial value is shown in table 3. It can be seen that, starting from 4 precession turns, the result of classification does not change, and only one chaotic orbit according to Minkowski dimension is classified as regular by $J_\parallel$
(this orbit is the same regular orbit mis-classified by Minkowski dimension above). It can also be seen that there are 8 regular orbits classified by condition (2.4) as chaotic. All these are regular transient orbits (see figure 2c) constituting less than 3 % of all regular orbits.
Table 3. Comparison of classification by Minkowski dimension and by using $J_\parallel$
. No more changes appear for more than 4 precessions, up to at least 128 precessions.
2.3. Summary of classification
Up to now, we have introduced two disjoint pairs of orbits classes: regular/chaotic and ideal/non-ideal. Ideal orbits are a subset of regular orbits, but non-ideal are not necessarily chaotic. Besides that, a common classification splits the trapped particle population into transient/non-transient. As mentioned above, transient means that an orbit undergoes transitions of the second kind between trapping classes during its precession, thereby abruptly changing its $J_\parallel$
value (see § 2).
Note, that in the bounce-averaged approximation (Velasco et al. Reference Velasco, Calvo, Mulas, Sánchez, Parra and Cappa2021), only transient orbits appear as chaotic. In the full guiding-centre model, where the bounce phase is retained, transient orbits cannot only be chaotic (as for the ones studied, e.g. in Kolesnichenko et al. Reference Kolesnichenko, Lutsenko and Tykhyy2022), but also regular. At the same time, ideal orbits are never transient. A graphical summary of this taxonomy with reference to examples is given in figure 6.
Figure 6. Taxonomy of orbits. Disjoint sets regular/chaotic, ideal/non-ideal and transient/non-transient are shown with their overlap and reference to example figures. The outermost regions are chaotic non-transient orbits.
From the analysis above we can conclude that in the considered test case in the QI configuration of Drevlak et al. (Reference Drevlak, Beidler, Geiger, Helander and Turkin2014), classification by $J_\parallel$
is more accurate than by Minkowski dimension with the current settings and much faster: it requires 4 precession turns i.e. approximately 100–150 footprints in contrast to near 10 000 footprints needed for the Minkowski dimension. The topological classifier with monotonicity reaches similar performance. A disadvantage of the $J_\parallel$
classifier is that the threshold has to be set by hand based on empirical data. Despite recognizing some non-ideal regular orbits, transient regular orbits are still classified as chaotic by the $J_\parallel$
criterion. The topological classifier only recognizes the subset of ideal orbits as regular, and thus requires more orbits to be traced to the end. In turn, when converged, the topological classifier requires no manual tuning.
It should be noted that a classification without errors by any method would require tracing orbits over infinite time. The main problem is presented by orbits close to periodic orbits (invariant axes, where the footprints are mapped onto themselves after a finite number of returns). According to the Poincaré–Birkhoff theorem (Lichtenberg & Lieberman Reference Lichtenberg and Lieberman1983) such orbits exist everywhere in the phase volume, including ergodic regions. Such regular orbits have very low rotational transform if they are close to an invariant axis. Thus a classifier can easily confuse them with ideal orbits, if the tracing time is too short. This behaviour is shown in figure 5. The closer an orbit is to an invariant axis, the longer it takes to correctly classify it. The same is true for chaotic orbits. In contrast to regular orbits, chaotic orbits depart from an invariant axis according to an exponential law (Lyapunov exponent). This departure can still be slow if they start very close to the axis. This can lead to a mis-classification as regular, but fortunately, the amount of these orbits is exponentially small.
3. Results and discussion
Figures 7–12 show the results from SIMPLE in the two new classifier modes as well as the collisional mode for five optimized stellarator reactor configurations (Subbotin et al. Reference Subbotin, Mikhailov, Shafranov, Isaev, Nührenberg, Nührenberg, Zille, Nemov, Kasilov and Kalyuzhnyj2006; Drevlak et al. Reference Drevlak, Beidler, Geiger, Helander and Turkin2014, Reference Drevlak, Beidler, Geiger, Helander, Henneberg, Nührenberg and Turkin2018; Henneberg et al. Reference Henneberg, Drevlak, Nührenberg, Beidler, Turkin, Loizu and Helander2019; Landreman & Paul Reference Landreman and Paul2022). For each configuration, classification runs are performed for 500 000 collisionless trapped alpha orbits started from banana tips ($v_\parallel =0$
) randomly distributed in plasma volume. Classification results are scored on the rectangular $50\times 300$
grid in starting banana tip radius $s$
and the normalized perpendicular invariant, $J_\perp =v_\perp ^2 B_{\rm min} / (v^2 B)$
which takes values ranging from zero for strongly passing orbits to $1$
for orbits deeply trapped at the global magnetic field minimum $B_{\rm min}$
. As a result, each grid cell contains $N_g$
good (ideal or regular, depending on classifier), $N_a$
average (non-ideal or chaotic) and $N_b$
bad (promptly lost) orbits. The predominant orbit class $C$
defined via these numbers as
is shown as a coloured pixel on this grid. Orbits are classified as promptly lost if they leave the confinement region before a criterion can be computed. If they reach this point, in the case of the $J_\parallel$
criterion, orbits are counted as regular, if the relative variation in $J_\parallel$
is below 1 %. Otherwise, they are counted as chaotic. For the topological classifier, orbits are counted as ideal if they fulfil the recurrence and monotonicity criteria, and non-ideal otherwise.
Figure 7. Quasi-isodynamic configuration of Subbotin et al. (Reference Subbotin, Mikhailov, Shafranov, Isaev, Nührenberg, Nührenberg, Zille, Nemov, Kasilov and Kalyuzhnyj2006). (a,b) show on an $(s,J_\perp )$
grid the predominant orbit class (3.1) in the trapped region bounded from above (deeply trapped $\boldsymbol {-}$
), and from below (trapped–passing boundary $\boldsymbol {--}$
). Thin vertical lines at $s=0.25$
mark starting flux surface for loss computation in (c,d). (a) shows $J_\parallel$
classification. Colour map for predominantly regular (bright), chaotic (medium) and promptly lost (dark) orbits. (b) Analogous to (a) for topological classification with ideal (bright), non-ideal (medium) and promptly lost (dark). (c) Left axis: $J_\perp$
for lost particles marked by dots and estimated loss distribution over time. Normalized average kinetic energy ($\boldsymbol {-\cdot }$
) of slowing-down particles. Right axis (red): collisional confined fraction of particles ($\boldsymbol {--}$
) and energy ($\boldsymbol {-}$
) over time. (d) Final lost energy fraction over $J_\perp$
for collisional (thick line) and collisionless (thin line) alpha particles. This configuration shows only minimal collisionless losses near $J_\perp = 0.55$
. This loss channel is marked by a white arrow in (b). Another loss channel near $J_\perp = 0.9$
is blocked for collisionless orbits, but not if collisions change $J_\perp$
over time. Thin dashed lines in (c,d) mark the minimum $J_\perp$
value from (a,b).
Figure 8. As in figure 7, but for the QI configuration of Drevlak et al. (Reference Drevlak, Beidler, Geiger, Helander and Turkin2014). Here, non-ideal orbits and losses extend over a region at the trapped–passing boundary near normalized $J_\perp = 0.7$
. The total collisional energy loss fraction slightly above $5\,\%$
matches the one in figure 4 of Drevlak et al. (Reference Drevlak, Beidler, Geiger, Helander and Turkin2014), where starting distribution and plasma profiles have been modelled in more detail.
Figure 9. As in figure 7, but for the QH configuration of Drevlak et al. (Reference Drevlak, Beidler, Geiger, Helander, Henneberg, Nührenberg and Turkin2018) with alpha losses via a region of non-ideal orbits around normalized $J_\perp =0.85$
.
Figure 10. As in figure 7, but for the QA configuration of Henneberg et al. (Reference Henneberg, Drevlak, Nührenberg, Beidler, Turkin, Loizu and Helander2019). Most trapped orbits are chaotic due to resonance overlap shown in figure 13(a), or promptly lost. Two regions with a larger fraction of regular but non-ideal orbits block some losses.
Figure 12. Quasi-axisymmetric configuration of Landreman & Paul (Reference Landreman and Paul2022) as in figure 7. The resonance structure matches figure 13(c) and the small deviation from quasi-symmetry prevents resonance overlap and therefore chaos. Prompt losses at $s\geq 0.5$
are from wide banana orbits.
To interpret the classification plots in the upper part of figures 7–12 it is important to recall from § 2 that, independently of the chosen thresholds, there exist regular non-ideal orbits. If the classifier makes no error, there exist, however, no chaotic ideal orbits. Completely regular regions and thus their subset of ideal regions guarantee collisionless confinement. Furthermore, collisionless orbits move on contours of constant $J_\perp$
. This means that whenever there exists a regular region with the same $J_\perp$
at a radius outside the initial position, this region will block the orbit from escaping. In the graphical representation interpreted as a projection in phase space this means that orbits may only evolve in the horizontal direction. Thus, regular/ideal regions plotted in bright yellow act as loss barriers. For the purpose of optimizing a stellarator configuration in terms of chaotic alpha losses, it should be sufficient to ‘pin down’ such barriers by maximizing the fraction of regular or better ideal orbits on a few flux surfaces, e.g. at $s=0.25$
and $s=0.5$
. The reason why ideal orbits are preferable to just regular ones will become clear below. From the computing time required for the present global analysis with the topological method, this should be feasible in a few minutes of computing time on a single cluster node and can be parallelized further easily.
3.1. The role of orbital resonances
Figure 11 is a special case to demonstrate the difference of lower-energy orbits of 35 keV to alpha orbits at 3.5 MeV in figure 10. At lower energy, the orbit width decreases, and the ideal orbit criterion converges towards the criterion of closed $J_\parallel$
contours (Gori et al. Reference Gori, Lotz and Nührenberg1996) in the limit of infinitesimal orbit width. The latter criterion is valid for thermal ions and usually fulfilled in a large fraction of phase space for devices close to omnigeneity (QI, QH and QA). In contrast, for energetic fusion alphas at 3.5 MeV, no predominantly ideal regions, and few regular regions, exist in this example. This behaviour is understood in terms of overlapping drift-orbit resonances computed in the code NEO-RT (Albert et al. Reference Albert, Heyn, Kapper, Kasilov, Kernbichler and Martitsch2016, Reference Albert, Rath, Buchholz, Kasilov and Kernbichler2022) and plotted in figure 13. An orbital resonance occurs whenever the canonical bounce frequency $\omega _b$
and toroidal precession frequency $\varOmega _t$
fulfil the condition
with integer harmonics $m_2$
and $m_3$
in canonical angle variables. This resonance is defined in the perfectly symmetric limit with three invariants of motion $\boldsymbol {J}$
– in this case quasi-axisymmetry. The special case $m_2=0$
where the orbit does not precess toroidally is called a superbanana resonance. Generally, in a configuration close to perfect symmetry, orbits perform nonlinear oscillations around such resonances that appear as island chains in Poincaré plots (non-ideal, see figure 2(a,b). The amplitude of these oscillations is known as the resonance width, and scales with the square root of the amplitude $H_{\boldsymbol {m}}$
of the Hamiltonian perturbation $\delta H$
in Fourier harmonics for canonical angles (Albert et al. Reference Albert, Heyn, Kapper, Kasilov, Kernbichler and Martitsch2016, Reference Albert, Rath, Buchholz, Kasilov and Kernbichler2022). In our case the perturbation scales linearly with the non-quasi-symmetric part $\delta B = B-B_0$
of the magnetic field module,
where quantities with subscript ‘0’ refer to their quasi-symmetric part.
Figure 13. Orbital resonance lines (Albert et al. Reference Albert, Rath, Buchholz, Kasilov and Kernbichler2022) for QA configurations from (a) figure 10, (b) figure 11 and (c) figure 12 in the thin orbit limit (Albert et al. Reference Albert, Heyn, Kapper, Kasilov, Kernbichler and Martitsch2016) of NEO-RT (precise values of $J_\perp$
and shape of resonances may differ from full orbits). Colours indicate different canonical bounce harmonics, and error bars show the linear resonance width.
Once the resonance width becomes of the order of the distance between resonances, chaotic phase-space regions appear in between. The degree of overlap is measured by the Chirikov criterion (Chirikov Reference Chirikov1979). Intact KAM surfaces in between island chains of resonances prevent crossing of chaotic orbits and are identified as footprints of ideal orbits.
For full alpha energy in figure 11, most orbital resonances overlap and produce global chaos. At reduced energy, only the superbanana resonance remains. Nonlinear oscillations around this resonance can be seen in the form of regular but non-ideal orbits in figure 11(b). Further away from the resonance line on the inboard side, one may even notice a region of prompt losses arising from nonlinear oscillations around the resonance over a radial region that crosses the plasma boundary $s=1$
. Such orbits are barely trapped in the resonance and perform wide oscillations over both, $J_\parallel$
and $s$
(figure 2b). Also a comparison of the second QA configuration in figure 12 with 13(c) shows a banded structure for isolated drift-orbit resonances, and a chaotic region from overlapping resonances near the trapped–passing boundary.
3.2. Classifiers as proxies for collisional losses
A comparison of $J_\parallel$
classifier with the topological classifier shows that the fact that regular orbits are a superset of ideal orbits is mostly reflected. Visually, this means, that the left plots show larger bright yellow regions than the right plots. An exception is seen, e.g. in figure 12 in the deeply trapped region. Here, the choice of a relative threshold of 1 % in $J_\parallel$
variations breaks down for very short banana orbits and leads to type I errors of this classifier that identify too many orbits as chaotic. This is proven by the topological classifier that identifies mostly ideal orbits in this region, and demonstrates that careful tuning of the $J_\parallel$
classifier is necessary.
The final judgment of the usefulness of the new classifiers and the kind of visualization discussed above is made based on collisionless and collisional direct loss computations in the lower panels of figures 7–12. Collisional computations are based on a relatively simple background model of constant temperature $T=10\,\mathrm {keV}$
for all species, and electron density of $n_e=10^{14}\,\mathrm {cm}^{-3}$
in a fusion plasma with a 50:50 D-T mixture. The left plots show the distribution of collisional alpha particle losses over time and initial $J_\perp$
. The right plots show the final lost energy fraction of collisional and collisionless computations. In all tested configurations, the loss channels of chaotic/non-ideal and promptly lost alpha orbits are reflected in the distribution of losses over $J_\perp$
.
4. Conclusion
From the investigations in this work, we can conclude that collisionless and collisional fusion alpha losses in stellarators are prevented by regular phase-space regions acting as radial barriers. In turn, transitions between trapping classes and overlap of drift-orbit resonances introduce chaotic regions that act as loss channels. In particular, overlap of resonances whose width scales with the square root of the perturbation amplitude leads to chaotic losses similar to those induced in tokamaks by toroidal field ripples. This analogy is seen in QA configurations and is expected to hold also in other quasi-symmetric as well as QI devices. Relevant regions in phase space can be efficiently charted by massive classification runs in SIMPLE after 1 % of the tracing time of a direct computation of loss fractions. In contrast to their superset of regular orbits that can also produce island chains in phase space, ideal orbits correspond to KAM surfaces that are more robust to perturbations. Ideal orbits are a natural generalization of closed $J_\parallel$
contours (Gori et al. Reference Gori, Lotz and Nührenberg1996). They can be identified at similar computational cost by the presented topological classifier without manual adjustments. We therefore propose to use the fraction of ideal orbits on a selection of flux surfaces as a metric for stellarator optimization.
In contrast to alpha particle confinement metrics based on the bounce-averaged motion (Gori et al. Reference Gori, Lotz and Nührenberg1996; Nemov et al. Reference Nemov, Kasilov, Kernbichler and Leitold2005; Velasco et al. Reference Velasco, Calvo, Mulas, Sánchez, Parra and Cappa2021; LeViness et al. Reference LeViness, Schmitt, Lazerson, Bader, Faber, Hammond and Gates2023), classifiers presented here are based on real orbits retaining their actual radial widths which can be compatible with plasma radius (see Albert et al. Reference Albert, Rath, Buchholz, Kasilov and Kernbichler2022; Landreman & Paul Reference Landreman and Paul2022) and, more importantly, bounce phases. Therefore, they naturally take into account drift-orbit resonances which may lead to chaotization of the orbits and respective additional losses. Such losses can be significant in stellarator configurations produced from real modular coils which result in short scale magnetic field modulations too small to create additional local wells (and thus significantly enhance bounce-averaged transport) but sufficient to generate resonant transport and respective losses. The later can be seen from the resonance condition (3.2) where the banana precession frequency $\varOmega _t$
being formally a higher-order quantity in Larmor radius than the bounce frequency $\omega _b$
can normally ($m_3 \sim m_2$
) match the resonance condition only for a small number of particles near class transition boundaries where $\omega _b$
drops. In case $m_3 \gg m_2$
realized for configurations where the precession direction is across the planes of modular coils (e.g. the case with toroidal field ripple in tokamaks) the resonance condition can be realized also for the bulk trapped particles. An exception are quasi-poloidally symmetric (Spong et al. Reference Spong, Hirshman, Berry, Lyon, Fowler, Strickler, Cole, Nelson, Williamson and Ware2001) or QI configurations where the precession direction is mostly within the modular coil planes. For the latter configurations even with real modular coils, as considered here (see also Sánchez et al. Reference Sánchez, Velasco, Calvo and Mulas2022), the effect of drift-orbit resonances is less significant. Nevertheless, the metrics proposed here can be useful for coil optimization also in those configuration types.
Acknowledgements
The authors thank P. Helander, M. Landreman, J. Loizu, M. Mikhailov and E. Paul for magnetic equilibria and useful discussions, together with the teams of the EUROfusion TSVV Task 12 on Stellarator Optimization and the Simons Collaboration on Hidden Symmetries and Fusion Energy.
Editor Per Helander thanks the referees for their advice in evaluating this article.
Funding
This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 – EUROfusion). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. We gratefully acknowledge support from NAWI Graz. The present contribution is supported by the Helmholtz Association of German Research Centers under the joint research school HIDSS-0006 ‘Munich School for Data Science - MUDS’.
Declaration of interests
The authors report no conflict of interest.
