Energetic particle optimization of quasi-axisymmetric stellarator equilibria

An important goal of stellarator optimization is to achieve good confinement of energetic particles such as, in the case of a reactor, alphas created by deuterium–tritium fusion. In this work, a fixed-boundary stellarator equilibrium was re-optimized for energetic particle confinement via a two-step process: first, by minimizing deviations from quasi-axisymmetry (QA) on a single flux surface near the mid-radius, and secondly by maintaining this improved QA while minimizing the analytical quantity ΓC , which represents the angle between magnetic flux surfaces and contours of J|| , the second adiabatic invariant. This was performed multiple times, resulting in a group of equilibria with significantly reduced energetic particle losses, as evaluated by Monte Carlo simulations of alpha particles in scaled-up versions of the equilibria. This is the first time that energetic particle losses in a QA stellarator have successfully been reduced by optimizing ΓC . The relationship between energetic particle losses and metrics such as QA error ( Eqa ) and ΓC in this set of equilibria were examined via statistical methods and a nearly linear relationship between volume-averaged ΓC and prompt particle losses was found.


Introduction
The confinement of energetic particles is important to the success of any fusion reactor. In contrast to axisymmetric devices such as tokamaks, stellarators do not inherently confine trapped particle orbits and thus tend to have worse neoclassical confinement, particularly in the low-collisionality regime [1][2][3]. As the collisionality of energetic particles is * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. significantly lower than that of thermal particles, stellarators must be optimized to achieve good confinement of energetic particles in particular [4]. Confinement of these particles is of special importance if stellarators are to be considered as potential fusion reactors, as alpha particles produced by deuteriumtritium fusion reactions, for example, must be well-confined in order for a plasma to reach ignition. In addition, both these fusion alphas as well as fast ions created by neutral beam injection and ion cyclotron resonance heating should be wellconfined in order to maintain high temperatures in the core and avoid damage to the walls and plasma-facing components (PFCs).
A stellarator equilibrium can be considered quasisymmetric (QS) if the magnitude of the magnetic field strength B = B(ψ, mθ − nNϕ), where ψ is the coordinate indicating the magnetic flux surface, θ and ϕ are the poloidal and toroidal Boozer angles, m and n are integers, and N is the number of field periods. In a stellarator, the ∇B and curvature drifts of trapped particles align with contours of B if they are very deeply trapped and/or if the stellarator equilibrium is QS. Thus, by improving quasi-symmetry, trapped particle drifts are aligned with flux surfaces and are better confined. In particular, in an equilibrium in which n = 0 and thus B is independent of the toroidal Boozer angle, known as a quasi-axisymmetric (QA) equilibrium [5], trapped particles drift toroidally rather than radially. These QA equilibria have trapped particle orbits approximating those of tokamaks.
It has been shown that QA can only be exactly achieved on, at most, a single flux surface [6,7], though it is possible to decrease the non-axisymmetric modes to a level lower than that of the geomagnetic field throughout a plasma volume, thus achieving very low energetic particle losses, at least in an equilibrium with zero plasma pressure [8]. Via the same method, configurations with very low deviations from QA and nonzero plasma pressure have been found, but only for ⟨β⟩ (the volume-averaged ratio of plasma pressure to magnetic pressure) up to 2.5% [9]. In the case of equilibria with higher ⟨β⟩, previous research has found that, by minimizing the deviation from QA (E qa ) on a single flux surface between half-radius and the edge of the plasma, energetic particle confinement in a QA stellarator can be significantly improved [10,11]. This can be easier and less computationally expensive than optimizing for energetic particle confinement directly.
Another analytic quantity that can serve as a proxy for energetic particle confinement in stellarators is the Γ C metric [12]. While minimizing deviations from quasi-symmetry seeks to align contours of B with flux surfaces, minimizing Γ C seeks to align contours of the second adiabatic invariant J || =¸v || dl with flux surfaces. This integral is performed between bounce points along a magnetic field line, with v || referring to the component of velocity parallel to the magnetic field. As collisionless trapped particles tend to follow contours of J || , reducing Γ C causes particles to drift along flux surfaces rather than away from them. This method has been used in quasihelically-symmetric (QH) configurations, resulting in significant improvements in energetic particle confinement [13].
In this work, this method was successfully applied to a QA stellarator for the first time. QH and QA stellarators tend to have different loss mechanisms which are differently correlated with Γ C [14], and so it was not obvious that the success of Γ C optimization in QH equilibria would translate to QA. The optimization done in this paper followed a two-step process: first, the starting equilibrium, LI383 (a three-period QA stellarator equilibrium that formed the basis of NCSX [15]) scaled down to a volume-averaged magnetic field strength ⟨B⟩ of 0.5 T and with a ⟨β⟩ of 4.3%, was optimized using the STEL-LOPT code suite [16] to have minimal E qa on a given flux surface. Following this, the resulting equilibria were optimized to minimize Γ C on several different flux surfaces (between 1 and 4) while maintaining low E qa .
Next, in order to determine the success of optimization, all resulting equilibria were scaled up to match the volume V = 444 m 3 and volume-averaged magnetic field strength ⟨B⟩ = 5.86 T of ARIES-CS [17], and Monte-Carlo simulations of fusion alpha particles in these devices were performed using the code BEAMS3D [18,19] in order to assess losses. Simulations were performed both with and without collisions.
Fifteen (15) equilibria, including several that were not viable candidates for actual fusion devices due to absence of favorable stability properties, were analyzed with statistical methods to find which analytical metrics best predicted confinement. Twelve (12) different metrics were calculated for each equilibrium, consisting of E qa , Γ C , and ϵ 3/2 eff (a term in the neoclassical transport coefficient in the 1/ν transport regime [3]), each averaged over various volumes within the plasma. Pearson's correlation r [20] was calculated for each metric with collisionless particle losses after 1 ms and 200 ms, as well as for collisional energy losses, in order to determine which metrics are the most useful for predicting energetic particle confinement.
In general, we would expect these three metrics to have similar relationships to energetic particle confinement. A perfectly QA device would have both ϵ 3/2 eff and Γ C = 0, and so, if we could reduce E qa to zero, it would be unnecessary to also target Γ C or to consider ϵ 3/2 eff . However, in practice, achieving perfect quasi-symmetry is not possible, and approaching the problem by minimizing Γ C might prove more fruitful, as while E qa = 0 means that Γ C = 0, the converse is not true, meaning that there are more possible equilibria to be found with minimal Γ C .
As for ϵ 3/2 eff , which was not targeted in the optimization, some correlation to energetic particle losses is expected, but previous work has shown that improving quasi-symmetry and Γ C can sometimes degrade ϵ 3/2 eff while still improving energetic particle confinement [13]. For this reason, we might expect ϵ 3/2 eff to be a weaker predictor of losses than the two targeted metrics. In this work, we sought to empirically test this and to determine which quantity of the three had the strongest correlation.
The layout of the paper is as follows. In section 2, we describe the methods and results of the two-step optimization for E qa and Γ C . In section 3, we describe the energetic particle simulation methods and results. In section 4, we define analytic metrics which are then statistically correlated with confinement. Finally, section 5 will give conclusions and implications for future optimization efforts.

Optimization
STELLOPT is a code which optimizes ideal MHD equilibria by varying the boundary terms used as input to the MHD equilibrium solver VMEC [21]. Many different algorithms are available for optimization: a modified Levenberg-Marquardt algorithm (bounded or unbounded), genetic algorithm with differential evolution, particle swarm evolution, and over 30 more, available through the suite MANGO 5 . The optimizations for this work were performed using both the bounded and unbounded Levenberg-Marquardt algorithms with finite difference derivative evaluation.
STELLOPT attempts to minimize a cost function where f target i can be MHD equilibrium quantities such as major radius R 0 , β, etc or non-linear functions of the equilibrium. The σ i are set by the user and determine the target weight w i = 1 σi , where higher weights mean more emphasis is placed on that target.
In fixed-boundary mode, which was used for all optimizations in this paper, the stellarator symmetric [22] VMEC boundary is represented by two sets of Fourier coefficients R m,n and Z m,n , defining the boundary as: in which m and n are the poloidal and toroidal wave numbers, θ and ϕ are the poloidal and toroidal angles, and N is the number of field periods. STELLOPT allows the use of different representations of the VMEC boundary. Our optimizations were done using the Garabedian representation [23], in which the boundary is converted into a complex number in which u and v are the normalized poloidal and toroidal angles: In addition, STELLOPT uses the normalized toroidal flux coordinate s: where Φ is the toroidal magnetic flux. In this work, the targets chosen for STELLOPT included the following equilibrium parameters: √ gdV/´√gdV, where p(s) is the plasma pressure, which was held constant during the optimizations, and µ 0 is the vacuum magnetic permeability Only equilibria in which these parameters varied less than ±5% from their original values were selected. Note that including the volume was redundant, as it is determined by R 0 and A. However, as it was used in the optimizations, it is included here for completeness. In addition, it should be noted that the toroidal current was held constant for all optimizations. Later in this work, the impact of the optimization on toroidal current profile is assessed and discussed.
In addition, two other parameters targeted in many of the optimizations were the second principal curvature κ 2 [24] and the magnetic well quantity. κ 2 is the minimum curvature of the plasma boundary, which was targeted in order to avoid an increase in the concavity of the plasma boundary.
The magnetic well target in STELLOPT is defined as as in [25]. (In this equation, ⟨. . .⟩ is used to denote a flux surface average rather than an average over the full plasma volume.) To have stability, we want this quantity to be positive on all flux surfaces. Note that this contains a term dependent on the plasma pressure p(V), and thus is not the same as the 'vacuum well' sometimes targeted in optimization. This version of W was selected as a target because the optimizations were done at finite β. It was found in experiments on LHD that despite the lack of a vacuum magnetic well, the pressure term stabilized the plasma [26]. Thus, this version of W, which includes the plasma pressure, gives a better picture of whether an equilibrium is stable at a given β. However, this does not mean the plasma will be stable at lower β, such as during startup, which was considered to be beyond the scope of this work. In some optimizations for this work, these two parameters were omitted. In others, they were included as optimization targets, either targeted to their original value in the LI383 equilibrium, or else given a floor via an altered version of χ 2 i . In the latter case, the typical cost function in equation (2) was replaced with a penalty function, making χ 2 i very large when the target dropped below the given target value and very small when the target was above it: where the constant a i and the weight 1 σi determine the steepness of the cutoff. Each of the three options were employed depending on which gave the best optimization results.

Quasi-axisymmetry optimization
For the first part of the optimization process, the above parameters (R 0 , A, V, β, W, and κ 2 ) were targeted while attempting to minimize QA error E qa . QA error in STELLOPT is defined as a measure of the prevalence of all non-axisymmetric modes of the magnetic field strength on a given flux surface. The QA error on that flux surface is defined as where the magnetic field strength B(s) on a surface is defined along the field line in straight-field line (or Boozer) coordinates [27] as a sum of Fourier coefficients B m,n : In the first part of our optimization, the goal was to minimize E qa as much as possible on a single flux surface. This was performed for a scan of five flux surfaces: s = 0.1, 0.2, 0.3, 0.4, and 0.5. It was not performed for any flux surfaces larger than s = 0.5 because as the optimization surface moved further outward, the shape of the plasma boundary became more complicated and thus less realizable by a finite coil set.
An additional optimization was performed, as it was discovered that most of the resulting equilibria had slightly negative W on the flux surface closest to the axis. The optimization of E qa on the s = 0.4 flux surface was repeated while increasing the steepness of the W cutoff (in other words, increasing a i and decreasing σ i in equation (10)) in order to avoid this issue, at the cost of some increase in E qa . This equilibrium was included with the others for analysis.
When this was finished, we were left with six new equilibria, each of which had E qa reduced as much as possible on a single flux surface without moving outside the boundaries for R 0 , A, V, β, κ 2 , or W, and without the plasma boundary becoming too complex.
In the rest of this text, these six equilibria will be referred to as QA1, QA2, QA3, QA4, QA4n, and QA5, with the number corresponding to 10s qa , where s qa is the surface at which E qa was minimized, and QA4n denoting the equilibrium found by increasing the weight of the cutoff for W. Shown in figure 1 are profiles of E qa (s) for all of these six equilibria, along with the starting equilibrium, LI383.
In addition, in figure 2, profiles of Γ C for each of these equilibria are shown. (Both E qa and Γ C are shown for s = 0.02, 0.04, . . . , 1.0.) We found that, simply by reducing E qa on a single flux surface, the Γ C metric had already been significantly reduced in all equilibria other than QA1, despite not being included as a target. (This is to be expected, given the relationship between E qa and Γ C discussed in the introduction: were E qa zero, Γ C would be as well.) The second step aimed to build upon this improvement by targeting this metric directly.

Γ C optimization
Γ C is an analytical metric defined in equation 61 of [12] that is used as a proxy for the energetic particle confinement of a given magnetic field configuration. In targeting Γ C in a QA device, the goal is the creation of toroidally closed contours of the second adiabatic invariant J || . Γ C depends on tan −1 vr v θ , where v r and v θ are the radial and poloidal particle drifts averaged over the bounce motion of the particles [12].
When Γ C is minimized, so is the ratio vr v θ , which means that confinement is improved either by minimizing v r , and thus the motion away from flux surfaces; by maximizing v θ , so that the poloidal drift of the particles leads them across the flux surface through areas of both positive and negative radial drift, averaging out the radial drift; or by a combination thereof. As ∂r , this process can also be described as one of aligning contours of J || with magnetic flux surfaces.
The Γ C metric was first used in the stellarator optimization code ROSE to improve energetic particle confinement in a QH stellarator. It was found that the largest improvements in energetic particle confinement came when simultaneously minimizing Γ C and deviations from quasi-helical symmetry [13].
In this work, our implementation of the Γ C metric in STEL-LOPT was utilized for the first time, and this was also the first time energetic particle confinement in a QA stellarator was successfully improved by optimizing Γ C . Rather than optimizing for Γ C and quasi-symmetry at the same time, we used a two-step process, beginning with the QA optimization. The implementation in STELLOPT was compared against the existing implementation in ROSE by calculating Γ C for the same VMEC input file, QA5, in both STELLOPT and ROSE, the results of which are shown in figure 2. The two codes agree relatively well, with the difference likely coming from the different integration methods used in calculating Γ C . The Γ C profiles calculated by STELLOPT are significantly less smooth compared to the E qa profiles. The reason for this is presently unknown, but it has been confirmed that it is not a matter of resolution, as the noise persists even at higher resolutions.
As will be discussed in section 3, the best performing QA-optimized equilibria were QA4 and QA5, and so these equilibria were chosen as the starting points for the Γ C optimization. An optimization from QA4n was done as well, in order to prioritize keeping W positive. Starting from these equilibria, we attempted to minimize Γ C on several flux surfaces simultaneously, while still staying within the limits for the other parameters (R 0 , A, V, β, W, and κ 2 ).
It should be noted that we also attempted to optimize the starting equilibria for E qa and Γ C simultaneously, as was done in [13]. This method produced less favorable results, in terms of reduction of the Γ C metric, than did the two-step process, and so the resulting equilibria were not included in this analysis. In the two-step process, E qa was included as a target in the second step in order to avoid minimizations of Γ C which might have increased E qa . However, weights were chosen such that the total χ 2 functional was to leading order determined by the minimization of Γ C rather than E qa or other targets. Because E qa was included as a target, it was sometimes further reduced in this second step of the optimization.
After optimization, Γ C was evaluated for each flux surface, giving a radial profile of this quantity for each equilibrium. While we also performed some optimizations in which Γ C was targeted on only a single flux surface, like E qa , we found that the best improvements in Γ C over the entire plasma were typically made when targeting three flux surfaces simultaneously, as was done in [13].
In the end, out of many optimization attempts, four configurations each were selected with optimized Γ C from starting equilibria in which E qa had been optimized at either s = 0.4 or s = 0.5. These equilibria will be referred to from now on as QA4-1, QA4-2, QA4-3, QA4-4, QA5-1, QA5-2, QA5-3, and QA5-4. The second number in each name refers only to the order in which the optimized equilibria were found. Note that QA4-1 was optimized starting from the equilibrium QA4n, while QA4- [2][3][4] were optimized from the starting point of QA4. Table 1 lists all of the flux surfaces on which Γ C was targeted for each of these eight optimized cases. One of the main differences in these cases was whether we attempted to minimize Γ C closer to the core or farther out; this was varied in order to determine which optimization strategy would lead to the biggest improvements in confinement. In addition to the flux surfaces of Γ C optimization, other differences between these optimizations include the degree to which the Γ C and E qa targets were weighted, as well as whether W and κ 2 were targeted and whether a traditional cost function or the penalty function shown in equation (10) was used. In addition, some optimizations were performed with the bounded Levenberg-Marquardt algorithm, while others used the unbounded version. The radial profiles of Γ C for all eight resulting configurations, as well as their beginning equilibria, are shown in figure 3. Almost all eight equilibria show improvements in Γ C over their starting equilibria. QA5-4 especially stands out as having the biggest improvements within the region s ∈ [0.4, 0.8]. The equilibrium QA5-3 was included despite having Γ C similar to or higher than QA5 on most flux surfaces in order to determine how much the increased Γ C in the core compared to QA5 would impact energetic particle confinement, as the other equilibria shared similarly low Γ C in the core.
Note that not all of these equilibria would be usable in a plasma device: some have negative W on some flux surfaces, while others have outer flux surface shapes that would not easily be created with coils. In figure 4, flux surface shapes at three toroidal angles (ζ = Nϕ = [0, π/2, π]) are compared between the original equilibrium LI383 and two wellperforming optimized equilibria, QA4-2 and QA5-4. Both optimized equilibria have more areas of concavity on their outermost flux surfaces, adding to the difficulty of creating these field configurations with finite coil sets.
Prior research has found that there are sometimes equilibria with significantly improved Γ C and/or quasi-symmetry (and thus better energetic particle confinement), but that are not usable due to lack of a magnetic well [8,28]. In order to make comparisons to those configurations, and to find how much improvement in energetic particle confinement was possible, we allowed for equilibria with negative W on some flux surfaces. In figure 5, it can be seen that, for example, in the equilibria QA4-2 and QA5-4, W is negative at the outermost flux surfaces.
Rotational transform ι was not targeted in the optimization, but rather was allowed to evolve freely. Because of this, many of the optimized equilibria have ι profiles differing significantly from that of LI383, with higher ι and lower shear (the  first radial derivative of ι). To consider how the changes in ι might affect the new equilibria, we consider which rational values ι passes through.
As all equilibria considered in this paper have three (3) field-periods, the values of ι at which islands are likely to form are those where ι = n m , where n is a multiple of 3. Islands are especially likely when the ι profile has low shear as it passes through a rational surface. In addition, because modes fall off as ( r a ) m , they will be smaller for high-m rational values closer to the core.
One incidental benefit of the optimization is that 12 of 14 optimized equilibria increased ι enough to avoid the ι = 3 6 rational surface. On the other hand, they now have lower shear as they pass through the ι = 3 5 surface, but they also cross this rational surface closer to the core of the plasma. Some equilibria now have ι values approaching 3 4 as well, but none cross it. In figure 6, ι profiles for LI383, QA4-2, and QA5-4 are compared. It can be seen that QA4-2 in particular avoids the 3 6 rational surface, though it comes close to 3 4 .

Energetic particle simulations
In this section, we compare simulations of energetic particle losses for fifteen (15) equilibria: the original equilibrium, LI383; the six QA-optimized configurations, QA1, QA2, QA3, QA4, QA4n, and QA5; and the eight equilibrium optimized first for QA and second for low Γ C , QA4-1, QA4-2, QA4-3, QA4-4, QA5-1, QA5-2, QA5-3, and QA5-4. The goal of this was to measure the losses for a large number of similar equilibria, all beginning from LI383, so that in the following section, we can then quantify the correlation between energetic particle losses and analytic properties such as E qa and Γ C . In order to assess the losses of fusion alphas, each equilibrium was scaled up to match the volume (444 m 3 ) and volumeaveraged magnetic field (5.86 T) of the ARIES-CS reactor, as was done for the comparisons in [28]. Monte Carlo simulations of fusion alphas were performed with the code BEAMS3D, which models the trajectories of charged particle gyro-centers in a 3D equilibrium field calculated using VMEC. BEAMS3D  can be used to model neutral beam deposition as well as collisionless trajectories and collisional slowing-down of energetic particles. The trajectories of the charged particle gyro-centers are followed using the drift kinetic equations: is the magnetic moment, and v || = d ⃗ R dt ·b( ⃗ R) is the component of the velocity parallel to the magnetic field [18]. The effects of collisions and pitch angle scattering are included between integration time steps [29].
BEAMS3D's deposition model has been validated [30] and a preliminary validation and benchmarking of the slowing down model has been done [29]. The simulations for this work were done using a new feature created for BEAMS3D that models fusion births and uses the birth rates to initialize particles for simulation [19].
The density and temperature profiles are mapped to the background computational grid in (R, Z, ϕ) using the VMEC equilibria. The reaction rate is then calculated on this grid. When particles are initialized, they are distributed equally on the grid with a randomized pitch angle (that is, the angle γ = cos −1 (v || /v) of the particle's velocity relative to the magnetic field). The weighting of each particle is then determined by the fusion birth rate at its initial location divided by the volume of its cell and the number of particles within. In this way, particles from the edge region, in which the birth rate is several orders of magnitude lower, can still be sampled [19].
This feature can simulate not only alpha particles born of DT fusion reactions, but also T, H, and He 3 born of DD fusion. However, the birth rates for the latter three particles are several orders of magnitude lower than those of fusion alpha particles, and thus in this work, we only considered the birth of fusion alphas.
All simulations were done with the same temperature and density profiles: defined on the grid s = 0, 0.2, . . . , 1.0. These values of temperature and density were chosen for easy comparison to simulations in [28]. As the plasma composition was set to be a 50-50 mix of deuterium and tritium, n D and n T were each half of the value of n e . Since impurities were not included in the simulations, the Z eff used was 1. The fusion rate calculated by BEAMS3D for the simulations can be seen in figure 7.
As the temperature and density profiles were chosen for comparison to prior research, they are not consistent with the pressure profile used to define the equilibria. In order to assess the impact on the configurations, the drift-kinetic solver the Stellarator Fokker-Planck Iterative Neoclassical Conservative Solver (SFINCS) [31] was used to determine self-consistent bootstrap current densities with these new profiles for the scaled-up versions of LI383, QA4-2, and QA5-4, depicted, along with the original current profiles from the initial scaledup VMEC equilibria, in figure 8. It can be seen that the higher pressure profile used results in an increase in bootstrap current magnitude, as expected.
Also of interest is the fact that the current density profiles calculated by SFINCS are similar for the three configurations, which shows that the optimization did not change the bootstrap current by a large extent. Defining the volume-averaged percent difference in ⟨J · B⟩ (averaged over each flux surface) between the optimized and original equilibria as where i is the equilibrium and N s is the number of surfaces included in the average, we end up with a difference of 6.92% for QA4-2 and 5.90% for QA5-4. Thus, the choice to perform the optimization at finite β with a fixed pressure profile and bootstrap current is justified. This is important because the bootstrap current is expensive to calculate, and in an optimization with many steps, it is preferable to either assume it remains constant, make use of a proxy rather than a full drift-kinetic solution, or only recalculate the bootstrap current after multiple optimization steps, as was done here. It has recently been determined that proxy formulae for bootstrap current used in tokamaks can also be accurately applied to QS stellarators [9], and one promising approach might be to make use of these proxy formulae followed by a full solution after optimization has been performed. Our results suggest that such a multi-fidelity approach is warranted in stellarator optimization at finite β.
Two simulations were run for each equilibrium: one including the effects of collisions, and one without. Approximately 128 000 3.5 MeV alpha particles were initialized for each simulation, with weights calculated via the fusion birth function, and followed until they either thermalized (for the collisional simulations), orbited in the plasma for 200 ms (for the collisionless simulations), or hit the last closed flux surface (LCFS), at which point they were considered lost.
The loss fraction of particles-and, for the collisional simulations, energy-as a function of time were calculated as such: where i indexes over all N particles, j indexes only over the n lost (t) particles lost at time t, w i is the weight of the ith particle, E j is the energy of each particle j at the moment it is lost, and E 0 = 3.5 MeV is the birth energy of each particle. Note that in the collisionless case, since the particles do not lose energy to collisions, E lost (t) = N lost (t).
In this paper, N lost will be used to refer to the particle losses at the end time of the collisionless simulations, and E lost will refer to the energy losses at the end time of the collisional solutions. In the former case, the simulations ended after 200 ms, while in the latter case, they ended after 120 ms, after which time all of the particles had been thermalized or lost for all equilibria.
In figure 9, results of these simulations are shown, with N lost from the collisionless simulations and E lost from the collisional simulations being plotted for each of the 15 equilibria. Of the cases which were only optimized for QA, we see that QA4 and QA5 are the best-performing equilibria. QA4 has the lowest N lost in the collisionless simulations (18.2%), but when collisions are included, QA5 has the lowest E lost (11.6%, as compared to 13.1% for QA4). This replicates the results from [11], which found the lowest losses when minimizing E qa at s = 0.4 or 0.5.
Note that QA4n, which was an attempt to improve the magnetic well of QA4 while maintaining good confinement, has significantly higher losses, particularly in the collisionless case (N lost = 23.0%). Though we see in figure 1 that the two equilibria have similar levels of E qa , QA4n has higher Γ C than QA4, as can be seen in figure 2, and this is likely the cause of the difference in confinement.
Next, looking at all 15 equilibria, it can be seen that the equilibria with the lowest N lost are QA4-2 and QA4-3, being the only equilibria with N lost < 10% (9.8% and 9.3%, respectively), while QA5-4 has higher N lost (16.4%) but comparably low E lost (7.0% for QA5-4, vs. 6.6% and 7.4% for QA4-2 and QA4-3). In fact, QA5 and the equilibria optimized from it, excluding QA5-3, have the lowest ratios of E lost : N lost , ranging from 0.42 to 0.61. This is because these equilibria lose particles more slowly, giving collisions time to free trapped  particles. In fact, other than QA1, every optimized case loses particles more slowly than the starting equilibrium LI383 and have a E lost : N lost ratio less than 1, showing that stochastic losses have been significantly decreased by the optimization.
In figure 10, N lost (t) and E lost (t) are plotted for LI383 and two of the best-performing cases, QA4-2 and QA5-4. It can be seen that, for times under 1 ms, N lost (t) in QA5-4 is less than or equal to N lost (t) in QA4-2, thus demonstrating the slower loss of particles that allows QA5-4 to perform better in the collisional case.
In addition, to explore the difference between the best and worst performing equilibria, a pair of particles were chosen from the collisionless simulations for LI383 and QA5-4. The particles were born at nearly identical locations (in (s, u, ϕ) coordinates, where u is the normalized poloidal angle in VMEC) and birth pitch angle. However, Particle 1 was quickly lost from LI383, while Particle 2 was confined in QA5-4. The orbit traces for Particle 1 and Particle 2 are plotted in (r/a = √ s,u) coordinates in figure 11. As mentioned previously, the goal of minimizing Γ C is to minimize the ratio vr v θ , where v r and v θ are the radial and poloidal particle drifts averaged over the bounce motion of the particles. It can be seen that this optimization was successful: in LI383, the banana orbit of the particle drifts radially, first inward and then outward, before reaching a location at which it becomes ripple trapped and drifts out of the plasma. However, in QA5-4, which was optimized for low Γ C , the banana orbit drifts in the poloidal direction and the particle remains well-confined.

Loss metrics
Once the energetic particle simulations were performed for all 15 equilibria, a statistical analysis was carried out to find the metrics which best correlated with energetic particle losses. First, we defined twelve analytical metrics which were expected to have some correlation to particle or energy losses, consisting of three different quantities each averaged over four different regions in the plasma.
Besides E qa and Γ C , which were targeted in the optimization, we also included the quantity ϵ 3/2 eff , which is a term in the coefficient of neoclassical transport in the 1/ν transport regime. For every equilibrium, each of these three quantities were calculated on 50 flux surfaces in the plasma, at s = 0.02, 0.04, . . . , 1.0. An averaged sum of each of these quantities was taken over four different regions: In order to evaluate which metrics were most strongly correlated to alpha particle losses, we calculated Pearson's correlation [20], r, between these twelve metrics and three different loss measurements: N lost , E lost , and N prompt = N lost (.001 s), or collisionless particle losses after 1 ms, also known as 'prompt losses'. Note that this is not the same as prompt losses in tokamaks, which refers to first orbit losses. Instead, N prompt measures deeply trapped particles which are very quickly lost, as in [32].
Pearson's correlation was calculated as: where the sample covariance is: and the mean X and standard deviation σ X are defined as: The value of r ranges from −1 to 1, with −1 and 1 representing perfectly linear relations (negative and positive, respectively) and 0 representing no correlation. In general, |r| > 0.5 is considered a large effect size. Note that r is unrelated to slope; it describes only how close the points are to falling on a line.
Here, X is one of the twelve analytical metrics, Y is one of the three loss metrics, the sample size, n, is 15, and i refers to each individual equilibrium. In order to get a larger sample, we included several cases which had been optimized, but were not usable due to issues such as negative W or plasma boundaries which are not realizable by a feasible set of field coils. All of the equilibria included were optimized from LI383 and maintained nearly the same aspect ratio, major radius, and β.
A previous study compared two similar metrics (both Γ C and deviation from quasi-symmetry evaluated at s = 0.6) with collisional energy losses, but the equilibria compared were very different from each other, including QA, QH, and non-QS equilibria [28]. This is the first time that these metrics have been compared to losses for a large group of very similar equilibria. Figure 12 plots all values of r for these metrics. The strongest effect sizes are seen for N prompt and E lost , rather than N lost . Looking at N lost (t) for various configurations, such as in figure 10, this seems to be because improving metrics such as  E qa and Γ C impacts not only the total losses N lost , but also how quickly those particles are lost. Thus, N prompt decreases as the metrics are improved, even if less improvement is seen in N lost . Similarly, E lost is impacted by changes in N lost as well as the time in which those particles are lost, as slower losses mean more time for collisions to have an effect, and so it can be a more useful gauge of improvements than N lost alone.
The strongest correlation found was between ⟨Γ C ⟩ and N prompt . The r calculated for these parameters was 0.995, a nearly perfectly linear relationship. In figure 13, these two quantities are plotted for all 15 equilibria, showing the strength of the relationship. ⟨Γ C ⟩ also had the strongest correlation with E lost of all 12 metrics (0.915), and the second strongest correlation with N lost , after ⟨Γ C ⟩ mid (0.500 vs. 0.568).
Note that ϵ 3/2 eff , particularly averaged over the entire volume and the edge, was also a strong predictor of losses despite not being targeted in the optimization. However, there was the confounding factor that the equilibria in this group with low E qa and Γ C also tended to have low ϵ 3/2 eff ; that is, it is unclear whether improving ϵ 3/2 eff alone, without corresponding improvements of E qa and Γ C , would have the same effect.
In addition, it seems that there is an inverse correlation between ⟨ϵ 3/2 eff ⟩ core and all measures of losses; however, as |r| < 0.2 here, this is not a strong relationship and may not hold true outside of this small sample of equilibria. The inability of ⟨ϵ 3/2 eff ⟩ core to predict energetic particle losses is in line with previous findings: it has been shown that configurations with improved Γ C and QS can have degraded ϵ 3/2 eff in the core, and that these equilibria tend to confine energetic particles much better than would be predicted from ϵ 3/2 eff alone [13].

Conclusions
In conclusion, it was possible to use a two-step optimization process, first minimizing deviation from QA (E qa ) on a single flux surface near mid-radius, and secondly minimizing the Γ C quantity on between 1 and 4 different flux surfaces, to obtain stellarator equilibria with lower energetic particle and energy losses. The lowest losses were found in equilibria in which E qa was minimized at s = 0.4 or 0.5, and then Γ C was minimized at s = [0.2, 0.4, 0.6] or s = [0.4, 0.6, 0.8]. Two configurations (QA4-2 and QA4-3) with N lost (collisionless particle losses after 200 ms) under 10% were found, and three (QA4-2, QA4-3, and QA5-4) were found to have E lost (collisional energy losses) between 6.6%-7.4%.
In contrast, the optimizations in this work were all performed at ⟨β⟩ = 4.3%.
In addition, in this work, the drift-kinetic solver SFINCS was used to find self-consistent bootstrap current profiles for the scaled up versions of the LI383, QA4-2, and QA5-4 equilibria, and the good agreement between the three equilibria suggests that stellarator optimization at finite β does not have a large impact on bootstrap current, thus justifying use of a multi-fidelity approach, such as keeping bootstrap current constant throughout the optimization and only performing the expensive full calculation after optimization, as was done here, or perhaps using proxies such as analytic formulae for bootstrap current used in tokamaks, which have been shown to be applicable to QS stellarators [9].
Would alpha particle losses of 7% be sufficiently low for a reactor? E lost also corresponds closely to the percentage of alpha power lost to the walls as calculated by BEAMS3D. Taking QA5-4 as an example, 638 MW of alpha power is born in the plasma. 6.7% (42 MW) is lost from the LCFS while the other 93.3% (595 MW) is absorbed by the plasma. Calculating the stored energy of the plasma from the profiles used in the simulations, we get 564 MJ, and dividing this by the absorbed power gives a required thermal energy confinement time (τ e ) of 0.95 s. This is similar to that required for ARIES-CS [17] and lower than values already achieved on the Joint European Torus [33], and so this can be said to be an achievable value of τ e . By this standard, we have achieved good energetic particle confinement for the purposes of reaching a burning plasma state.
An additional concern in energetic particle confinement is that of damage to the walls and PFCs due to impinging energetic particles. As mentioned above, in the collisional simulation for QA5-4, 42 MW were lost out of a LCFS surface area of 724 m 2 . However, these losses are not evenly distributed, and the risk of damage from lost energetic particles depends not only on the total loss but on the spatial distribution of these losses. Determining the loads to the wall and PFCs would require full-orbit simulations performed with a detailed wall model.
It should be noted as well that these equilibria are not reactor designs. One weak point of this work is the fact that the temperature and density profiles were chosen for comparison to prior research, rather than for self-consistency with the pressure profile and β. The 11.5 keV temperature value at the core is a pessimistic assumption; with this level of heating, the actual temperature might be twice that, and the density significantly lower. Despite these shortcomings, the results show the success of our optimization methods.
Compared to E qa , which was specifically targeted in the optimizations, and ϵ 3/2 eff , which was not targeted but is known to influence confinement, Γ C was the strongest predictor of losses. In particular, there was a nearly perfectly linear relationship between ⟨Γ C ⟩ (that is, Γ C averaged over the entire volume) and N prompt (prompt losses). There were smaller, but still significant, correlations between Γ C and E lost and N lost .
It is unsurprising that minimizing Γ C would strongly improve prompt losses, as it has been found that, for stellarators approaching quasi-symmetry, the most prevalent mechanisms for losses under 1 ms (drift-convective banana losses and ripple trapping losses) are driven by large positive values of Γ C . The significant but weaker correlation between Γ C and longer-timescale losses of both collisionless particles and energy is also in line with the properties of LI383 in particular, for which it has recently been shown that particles lost on longer timescales spend less time in drift-convective or ripple orbits, but most often are lost by eventually drifting into a ripple orbit [14].
This suggests that optimizing Γ C , or similar quantities such as Γ α or Γ δ , defined in [34], is a promising avenue to obtaining stellarator equilibria with good energetic particle confinement.