Simulating fusion alpha heating in a stellarator reactor

Gyrocenter following simulations of fusion born alpha particles in a stellarator reactor are preformed using the BEAMS3D code. The Wendelstein 7-X high mirror configuration is scaled in geometry and magnetic field to reactor relevant parameters. A 2×1020 m−3 density plasma with 20 keV core temperatures is assumed and fusion birth rates calculated for various fusion products assuming a 50/50 deuterium-tritium mixture. It is found that energetic He4 ions comprise the vast majority of the energetic particle inventory. Slowing down simulations of the He4 population suggest plasma heating consistent with scaled energy confinement times for a stellarator reactor. Losses for this configuration appear large suggesting optimization beyond the scope of the W7-X device is key to a future fusion reactor. These first simulations are designed to demonstrate the capability of the BEAMS3D code to provide fusion alpha birth and heating profiles for stellarator reactor designs.


Introduction
The heating and confinement of the fusion products in a stellarator reactor is explored using the BEAMS3D stellarator energetic particle code [1]. Stellarator experiments rely on external heating methods which are used to mimic the heating conditions of higher magnetic field, larger, energy producing reactor environments. However, unlike the reactor, these systems do not necessarily rely on the slowing down of energetic fast ions for heating of the plasma. And, while these systems are capable of producing energetic particles, it is not clear that their * 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. birth profiles are similar to those of fusion alphas in the reactor scale devices. In this work, simulations of alpha birth profiles are performed for a reactor scaled device and compared to neutral beam injection in the Wendelstein 7-X stellarator (W7-X).
Modern fusion devices are routinely heated by a combination of ohmic, neutral beam, and radio frequency heating methods, while fusion reactors are envisioned to have their plasmas predominantly heated by the slowing down of fusion born alpha particles. It is natural to ask how well do such auxiliary heating methods mimic that of alpha heating. To understand this question one must understand that the alphas will collide with the thermal background plasma imparting energy to both electrons and ions. As the 3.5 MeV birth energy is one to two orders of magnitude larger than the background plasma thermal energy, these particles will predominantly slow down on electrons. Thus it is correctly argued that direct heating of electrons via such auxiliary heating systems is the most relevant heating scenario to mimic a reactor scenario. In a stellarator (and in perturbed tokamaks), not all energetic particles will be confined and, as such, modeling of the particle orbits with slowing down and pitch angle scattering is necessary to estimate the fusion alpha heating.
In this work, we document the implementation of fusion alpha birth profiles in a gyrocenter following code capable of treating the three dimensional nature of stellarator magnetic fields. In section 2, the implementation is described along with our choice of scaled up W7-X equilibrium. In section 3, the simulations are presented with focus on birth profiles, particle orbits, slowing down, heating, and power loss through the equilibrium boundary. Finally, in section 4 we conclude the manuscript with a comparison between the fusion heating profiles and those obtained through auxiliary heating on W7-X.

Methods
The high mirror magnetic configuration of the W7-X stellarator is scaled up to reactor relevant parameters and fusion alpha confinement is evaluated with the BEAMS3D code. This magnetic configuration is considered as it is the quasiomnigenous configuration of W7-X [2]. As such it is the realization of the optimization criterion for W7-X and should exhibit favorable energetic particle confinement. The BEAMS3D code is used as it is interfaced to the VMEC equilibria [3] and previously validated against measurements in W7-X [4,5]. A newly developed capability in the code is used to model fusion product births. Slowing down simulations are preformed with the code to the plasma edge as defined by the VMEC code. Regions outside the equilibrium are not considered in this work as here the gyrocenter approximation does not produce a reasonable estimate of wall loads.
An enlarged version of the W7-X high mirror magnetic configuration is used in this work (figure 1). The basis equilibrium was constructed by fitting a VMEC boundary to the flux surface inside the n/m = 5/5 island chain as computed in vacuum. Here 128 toroidal flux grid points, 13 poloidal modes (m = [0, 12]), and 21 toroidal modes (n = [−10, 10]) are considered in the calculation. This equilibrium was then enlarged by multiplying the boundary harmonics by a constant, thereby preserving the aspect ratio while enlarging the plasma. This resulted in the plasma volume going from ∼26 to ∼1700 m 3 and the major radius going from ∼5.5 to ∼22 m. The enclosed toroidal flux was changed from ∼1.7 to ∼64 Wb taking the W7-X axis magnetic field from 2.5 to 5 T. The plasma profiles were assumed and equilibrium pressure chosen to be consistent, resulting in a plasma β of 3%. The plasma current profile was chosen to be zero for this work as it is expected to have little impact on the equilibrium used in this study. As a test the spectrum of | ⃗ B| in Boozer coordinates was checked and the new equilibrium had a similar structure to that of the W7-X high mirror case (other quantities like rotational transform were also preserved).
In order to model the confinement of fusion products arising from nuclear reactions a fusion birth model was implemented in the BEAMS3D code. The model is based upon the fusion cross sections presented in the paper by Bosch and Hale [6], which are a fit to experimental data. The fit takes the form of the equation: where C k are coefficients for different reactions, T i is the ion temperature (in keV), B g is the Gamov constant, m r c 2 is the reduced mass of the particles (in keV), n D is the deuterium (D) isotope density, and n T is the tritium (T) isotope density. In this work it was assumed that n D = n T = 0.5n e , which was done for simplicity and could be relaxed given profile modeling. Thus the background plasma was assumed to be a 50/50 mix of deuterium and tritium. In BEAMS3D, the equilibrium quantities are mapped to a cylindrical (R, ϕ, Z) grid and the ⟨σV⟩ particle birth rate is also mapped to that grid. The follow fusion reactions are considered in this work: and D + D → He 3 (0.82 MeV) .
In this work we only consider the reactions between thermal species and neglect those created through the thermal deuterium and fast tritium reactions. In Monte-Carlo codes such as BEAMS3D, markers are initialized to represent ensembles of energetic particles and followed. In such simulations the 'weight' of a marker acts to define the flux of particles each marker represents. For neutral beam simulations, fast ions are born by following markers along the neutral beam path and using the background charge exchange collisional cross sections as a way to determine the starting location for the fast ion. In that model the weight represents the fraction of particles injected by the neutral beam at a given energy. For nuclear reactions, the whole volume of the plasma gives birth to energetic particles. In this work a marker is initialized in every cell of the cylindrical grid which falls inside the equilibrium domain. It is then given a 'weight' according to that grid cell's volume and ⟨σv⟩. If more markers are requested than needed to evenly fill this space, the additional markers are seeded randomly. Cells with more than one marker, have the enclosed marker's 'weights'  Note that it is assumed that Z eff = 1 and Er = 0 in this work.

Quantity
Value Units 135 168 none suitably reduced. Initial pitch angles are randomly chosen for all particles. Table 1 highlights the simulation grid parameters. Initial attempts to generate markers based on reaction coefficients proved problematic. First, one would like every position which can birth to particles to have at least one marker. This places a minimum threshold on the number of markers to be the number of computational 'cells' in a given simulation. For a stellarator this can easily exceed 100 000 markers. Second, regions of high reaction rate may have cross sections which are order of magnitude larger than the lowest reaction rates. Thus many millions of markers would be needed to properly fill this space. One could place a lower threshold on the reaction rates considered but this quickly causes markers to not be born anywhere in the 'edge' region of the equilibrium. One can check the validity of such approximations by calculating the amount of power born in the markers against the amount of power which would be born in the total plasma (using integrals over the plasma volume and radial profiles). This comparison only converged for very large numbers of particles. The practice of a uniform grid of markers with varying weights proved to me the most computationally efficient method for representing the power born in the plasma.
It should also be noted that fusion births are isotropic in pitch angle space. For this work the pitch angle of each marker is randomly chosen. As the number of markers increases in a given 'cell' this space will become better filled. A detailed investigation of such phenomena is left to future work where computational efficiency of the code has been improved.
For simplicity two simulations are considered in this work. One in which all four fusion products are considered and one in which only the fast He 4 is considered. In both simulations slowing down and pitch angle scattering operators are employed, however the radial electric field is neglected (for lack of an appropriate model potential). In the first simulation, particles are followed for 100 time steps for 100 ns. Such simulations, allows the analysis of the relative importance of each species along with obtaining characteristic orbit information. The second simulation follows the fast He 4 only over a slowing down time (500 ms, as obtained from flux surface slowing down simulations). Justification for this down selection is made in the next section. Clearly the neglect of the fast tritium, fast protons, and fast He 3 species in the more detailed analysis is justified by this analysis. We point out that these are only the thermal births where we have assumed no beam-plasma interactions. Such an assumption is justified if one assumes no neutral beam injection will be present in the reactor. We also do not consider the so-called 'burn-up' of the fast tritium [7]. Such reactions would give rise to additional fast He 4 . Their neglect is justified given the ratio of birth rates. Figure 3 depicts the initial orbits for the different species over a time period of 50 µs. Passing orbits are generally well confined in toroidal systems so we will not focus on those here. Instead we have chosen to focus on two trapped orbits which are more stereotypical of W7-X fast particle orbits. The first  orbit considered is the poloidally precessing deeply trapped particle (in red, arctan

Results
. Particles such as these mirror inside a field period while precessing poloidally. Over time this poloidal precession moves the particle from bad to good curvature regions averaging out the radial drift of the particle. The second orbit is also a trapped orbit but closer to the trapped/passing boundary (in blue, arctan . This particle also bounced between mirror points but the mirror points are located in adjacent field periods. Such behavior is attributed to the fact that in the W7-X magnetic configuration the magnetic field strength is both a function of poloidal and toroidal angle. It should be noted that it is also possible for a poloidally precessing passing particle just outside the trapped-passing boundary, to become trapped for a time. This is again because of the nature of the three dimensional magnetic field. Such features highlight the more complex nature of orbits in stellarators. Figure 4 depicts plasma heating based on the assumption that the fast ions slow down on the flux surface on which they are born. Such an analysis neglects the effects of particle radial drifts and power which will be lost to the wall through imperfect confinement. However, it is useful for quickly checking the relative heating of the different fast ion species. Clearly the fast He 4 is the majority contributor to the heating, as such we will focus on this population only from this point (which also reduces the computational load of the simulations). From this analysis, an initial estimate of the time to slow down (τ )  of each species can be made: τ He 4 = 450 ms, τ T = 720 ms, τ p = 760 ms, and τ He 3 = 150 ms. This places an upper limit of how long to run a gyrocenter following simulation, which was set to 500 ms in practice.
The evolution of marker losses can be seen in figure 5 helping to show the interplay between particle losses and thermalization. The particle losses appear to begin around 100 µs into the simulation appearing to asymptote at around 1 ms. This is then followed by an increase in losses again around 5 ms into the simulation. Particle thermalization begins to rise around 10 ms into the simulation with all particles being either lost or thermalized by 460 ms into the simulation. Overall the energetic particle loss fraction is rather high. This a consequence of simply taking the W7-X high mirror configuration and scaling it up. Optimized Helias reactor configurations are the subject of ongoing works, which suggest such losses can be significantly reduced. Still this analysis suggests that collisionless simulations may result in unrealistically high loss fractions if followed longer than the plasma slowing down time.
When we compare the losses for this collisional run with a collisionless run we find two features. First the prompt losses (∼100 µs) are well resolved by both simulations. This population of particles is essentially behaving in a collisionless limit as is assumed by most works on stellarator fast ion optimization. The second feature is a strong disagreement on losses on longer timescales. The collisionless run sees a much smaller fraction of losses which has yet to asymptote at the 500 ms point. The enhanced losses in the collisional run are attributed to the fast ions slowing down and beginning to pitch angle scatter into the loss cone. This second point has strong ramifications for the assessment of optimized stellarator configurations.
In previous works, collisionless energetic particle confinement is used as a metric in the assessment of stellarator designs. Such simulations are usually run till the collisionless losses asymptote, usually in the ∼1 s range. However, this may not be a useful assessment as this is both longer than the stereotypical slowing down time and focuses the result on longer timescale orbits whose confinement is not well captured by collisionless orbit physics. At the very least such simulations should be viewed in terms of their collisional slowing down time rather than when collisionless orbits asymptote.
The heating profile arising from the slowing down simulations is plotted in figure 6. Here we can again see the electron heating dominating over the ion heating while the amount of heat deposited at all core radii appears to be reduced as compared to flux surface slowing down simulations. Integrals over these profiles indicate 60 MW of power going to the ions, 280 MW going to the electrons, and 40 MW being lost though the plasma boundary. Such analysis suggests that up to 25% of the fusion power is unconfined. One should remember that we have only scaled up the W7-X equilibrium for this work and as such it is not optimized for fast ion confinement. Other works are suggesting stellarators with very small energetic particle losses can be optimized [8][9][10][11].
It should be noted that here we have calculated the power lost by considering both the energy of markers as they hit the equilibrium boundary and their associated weights of the markers. If we were to only consider their weights then we'd find 34% of the particles are lost. Counting only the lost markers we'd find 36% lost. This highlights the need to perform such simulations taking care to account for the energy loss to the plasma on the way to the wall. Figure 5 clearly confirms this as only a small fraction of the particles are promptly lost at around 100 µs. As the weights of the particles account for births and weight the core markers over the edge markers, we can see that the majority of these losses are from the core.
Before continuing, we can estimate if this level of heating is sufficient to sustain the plasma. In this case our equilibrium stored energy is around ∼1 GJ which is equally distributed between the ions and electrons. Thus the total energy confinement time of the thermal plasma must be at least 2.6 s to sustain these plasma parameters with this level of heating. While such levels of confinement have yet to be achieved in smaller devices, they are not outside the scope of confinement times predicted when using scaling laws for stellarator reactors. Of course this discussion eludes the question of how said profiles are achieved and whether they are consistent with transport theory. Such analysis is left to future works. It is worth noting that earlier simulations with the 'so-called' narrow-mirror W7-X configuration (scaled as was done here), would only need confinement times of around ∼1 s which would significantly reduce the burden of finding enhanced confinement above scaling laws.
Examination of the lost particles suggest the presence of both prompt and longer timescale losses (figure 7). Here we define prompt as particles lost near their birth energy (3.5 MeV). Although prompt losses are clear in the plot, the majority of losses arise from particles which have significantly slowed down. Despite a uniform birth distribution in pitch angle space, lost particles have pitch angles inside |v ∥ /v total | = 0.5. This suggests that slowing down and pitch angle scattering processes are playing a large role in the losses. This is consistent with the picture of some passing particles slowing down, becoming pitch angle scattered into the loss cone, and finally being ripple lost. This is also consistent with the temporal evolution of losses. From 100 µs to 1 ms we see losses asymptote, these are the prompt losses, with losses growing from around 10 ms to the end of the simulation around 500 ms. Finally, we note that losses appear asymmetric in pitch angle space. The prompt and pitch angle scattered losses appear in the field aligned direction (v ∥ > 0), while a subtle band of losses at all energies appears in the anti-parallel direction (v ∥ < 0).
A final analysis of the heat flux through the plasma boundary (as defined by the equilibrium last closed flux surface) is shown in figure 8. Here we clearly see an up/down asymmetry is present in the losses along with a toroidal asymmetry, indicating that energetic particle orbits violate the flip mirror symmetry of the stellarator (stellarator symmetry). The highest heat loads appear on the upper outboard side of the bean shaped cross section and shifted slightly in the negative toroidal direction. While not directly indicative of wall loads, the loads are reaching above 10 MW m −2 . Similarly large heat loads are seen at the bottom inboard side of the triangular cross section. While only a small stripe is seen at a similar location on the top of the triangular cross section of the plasma. This result suggests that even relatively slow lost fast alphas can contribute to significant heat loads to the first wall. Modeling of actual wall loads is left to future work where full-orbit models and scrape-off layer effects can be assessed.

Conclusions
In this work the confinement, heating, and loss of fusion born fast alphas in a reactor sized stellarator equilibrium are investigated with the BEAMS3D code. The ability to model fusion products leverages the newly developed capability of the BEAMS3D code to generate fusion born energetic particles using realistic (if not optimistic) plasma profiles. Birth profiles show that the majority of heating, as long assumed, comes from the deuterium/tritium born fast helium-4 as opposed to fast protons, tritium, or helium-3. The comparison between the orbit following and simplified, slowing down on a flux surface, calculation shows that the inclusion of orbit dynamics reduces the overall deposited power to the ion and electron channel. However, losses are still unacceptably high (>10 MW m −2 ), a situation which can be rectified though further optimization of the magnetic configuration. Following particles to the equilibrium boundary, lost particles are found not to respect the flip-mirror symmetry of the stellarator magnetic fields. This implies that armoring of a reactor for lost fast ions may be possible, as only specific locations are loaded. Such a demonstration lays the groundwork for assessing future optimized stellarator reactor designs.
Before ending this manuscript it is worthwhile to place the alpha heating profiles of this scaled up W7-X in the context of the heating systems of W7-X. Here the quantity to compare is the power density as this is relatively device independent. In figure 9, we compare the alpha heating profiles to those achieved with the neutral beam heating. Four hydrogen sources with ∼1.8 MW of neutralized power are considered in this analysis. The ion heating appears consistent with power densities achievable through nuclear fusion, however a clear shortfall in electron heating is present. This is attributed to the 55 keV injection energy of the W7-X neutral beam system. The necessity of additional electron channel heating (in the form of electron cyclotron resonance heating) is clearly evident.
In this work, the effects of time varying fields, gryoorbit, charge-exchange, and scrape-off layer effects have been ignored. While it is possible that mode activity can perturb the magnetic field in stellarators, it is expected that such effects will be much smaller than those found in tokamaks. This is Figure 9. Comparison between alpha heating in the scaled up reactor to the NBI heating system (four sources). The ion heat deposition appears similar but a shortfall in electron heating is clearly present. due to suppression of Alfvénic activity [12,13], suppression of neoclassical tearing modes [14], and mitigation of sawteeth present in stellarators [15,16]. The high magnetic field (5 T) of this equilibrium justifies the use of a gyrocenter following code through a small gyro-radius. However, simulations extending to the first wall should incorporate full orbit following methods. The neglect of neutral physics in this work implies that the values from such simulations are upper limits. Work is underway to implement and validate change exchange models for fast ion confinement in stellarators. Future works will work to include edge effects and modeling of fast ion wall loads in optimized reactor scenarios.

Data availability statement
The data generated and/or analyzed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.