Study of Alfvén eigenmode stability in Quasi-Poloidal Stellarator (QPS) plasma using a Landau closure model

The aim of this study is to analyze the linear stability of Alfvén eigenmodes (AE) in the QPS device heated by a tangential neutral beam injector (NBI). The analysis is performed using the gyro-fluid code FAR3d, that solves the reduced MHD equations for the thermal plasma coupled with moments of the kinetic equation for the energetic particles (EP). The AE stability is calculated in several operational regimes of the tangential NBI: EP β between 0.001 and 0.1, EP energy between 12 and 180 keV and different radial locations of the beam. The analysis is performed for vacuum and finite β equilibria as well as QPS configurations with two and three periods. The EP β threshold in the vacuum case is 0.001 and the AE frequency is lower as the energy of the EP population decreases. Toroidal Alfvén eigenmodes with f = 80–120 kHz and elliptical AE between f = 120–350 kHz are triggered between the middle-outer plasma region ( r/a>0.5 ). The AE stability improves in the simulations with finite β equilibria and three period configurations with respect to the vacuum case with two periods because the continuum gaps are slender, leading to a higher threshold of the EP β, above 0.03 for the AEs triggered by the helical mode families. Helical effects are not strong enough to destabilize Helical Alfvén eigenmodes, the AEs with the largest growth rates are triggered by the n = 1 and n = 2 toroidal families.

The aim of the study is to analyze the Alfvén eigenmode (AE) stability of QPS plasma heated by a tangential neutral beam injector (NBI). The analysis assumes a tangential NBI identical to the tangential NBI in the LHD device [47] that injects EPs with a birth energy of 180 keV. The AE stability is calculated for QPS vacuum configurations with two and three periods as well as a QPS finite β scenario with two periods. The optimal NBI operational regime for improved AE stability is identified with respect to the NBI injection intensity (EP β) (considering the EP slowing down process) and the NBI deposition region. The stability of dominant and sub-dominant modes is calculated.
The analysis is performed using the gyro-fluid code FAR3d [50,[63][64][65]. FAR3d solves the reduced linear resistive MHD equations, adding the moment equations of the energetic ion density and parallel velocity reproducing the linear wave-particle resonance effects required for Landau damping/growth [66][67][68][69]. The set of equilibria used in the study is calculated by VMEC code [70]. A parametric analysis is performed with respect to the EP energy, EP β and EP radial density profile, including the effect of the helical couplings and finite Larmor radius (FLR) stabilization. This paper is organized as follows: the numerical scheme, equilibrium properties and simulation parameters are shown in section 2. The simulations results are described in section 3 for the two period vacuum, three period vacuum and two period finite β configurations. The conclusions and discussion are presented in section 4.

Numerical model
The FAR3d numerical model for the thermal plasma uses a set of reduced MHD equations retaining the toroidal angle dependency in a three-dimensional equilibrium [70,71]. The resonance of the EPs with the marginally stable shear Alfvén waves is introduced by moments of the gyro-kinetic equation truncated with a closure relation. The closure coefficients are optimized to bring in the effects of the parallel Landau resonance (Landau coefficients) [66] in the zeroth order (EP density n f ) and the first order (EP velocity parallel to the magnetic field lines v ∥,f ) moment equations. Landau closure coefficients used for the model calibration are calculated by matching the analytic TAE growth rate of the two-pole approximation of the fluid plasma response function with a kinetic response function for energetic particles obeying a Lorentzian distribution. The lowest order Lorentzian can be matched to a Maxwellian or to a slowing-down distribution by choosing an equivalent average energy (a Maxwellian for the FAR3d model). The FAR3d model includes a velocity averaged grad-B orbit drift propagator in each of the fast ion moment equations. This includes radial, poloidal, and toroidal grad-B drift components. Thus it effectively incorporates the average orbit deviation from magnetic surfaces by factoring these drifts into the fast ion moment dynamic evolution.
The simulations have been performed including the effect of the continuum and FLR EP stabilization effects [3]. The electron-ion Landau damping is not added because the effect is negligible due the low temperature of the plasma. EP FLR stabilization effects improves the stability of the modes reducing the available free energy required to destabilize the AEs. This effect leads to a decrease of the AE growth rate or the AE stabilization. The simulations in the present study only reproduces the destabilizing effect of passing EPs. The destabilizing effect of trapped EP is not included in the analysis. The pitch angle of the EP generated by a tangential NBI is small, thus the model is dedicated to studying the stability of AEs induced by passing EPs. Spong et al [67,68] include further details about the model equations and numerical scheme.

Equilibrium properties
A set of fixed boundary VMEC equilibria [60] are calculated for three QPS configurations with a hydrogen plasma: vacuum and finite β cases for a two magnetic field period configuration (called 2P Vacuum and 2P Finite β) as well as a vacuum case for a three magnetic field period configuration (called 3P Vacuum). The averaged inverse aspect ratio is ϵ = 0.257, the major radius R 0 = 0.9 m, the minor radius < a >= 0.35 and the magnetic field at the magnetic axis is 1T. The energy of the particles injected by the hydrogen NBI is a hypothetical T f = 180 keV [38] resulting in an averaged Maxwellian energy equal to the average energy of a slowing-down distribution of 90 keV.
The Larmor radius is calculated using the following expression: with m i is the ion mass, V th,f the average velocity of the EP distribution function, q the ion charge and B the magnetic field. The Larmor radius calculated for EP populations with energies between 30-180 keV ranges from 0.01 m to 0.05 m. The Larmor radius used in the simulations is 0.02 m. Figure 1 shows the main profiles of the thermal plasma. Thermal plasma temperature profile shown in panel (d) is the same for all the FAR3d simulations performed. It should be noted that the VMEC equilibria of the vacuum cases are calculated setting the thermal plasma pressure/temperature to zero. The AE stability of each model is mainly determined by the magnetic field configuration, reason why the same thermal plasma profiles are used in all the simulations for simplicity. The continuum gaps are calculated using the code Stellgap including the effect of the sound wave (using the temperature profile indicated in figure 1 and helical couplings for all the cases) [38].

Simulation parameters
The dynamic and equilibrium toroidal (n) and poloidal (m) modes included in the simulations are listed in table 1. The notation used to identify the modes is n/m. The linear simulations are limited to n = 1 to 5 toroidal modes. The mode selection includes all the resonant modes between the magnetic axis  [11,26] and the plasma periphery. The number of points of the radial grid is 1000. The dynamic variables must include both mode parities because the moments of the gyro-kinetic equation break the MHD symmetry. Eigenfunctions in FAR3d code are represented in terms of sine and cosine components, using real variables: The cosine and sine components of the eigenfunction are indicated by positive and negative toroidal mode numbers, respectively.
The EP populations leading to the strongest resonance are identified with respect to the ratio between the EP velocity (proportional to the square root of the EP energy) and the plasma Alfvén velocity. The thermal ion density and magnetic field intensity are fixed in the models, thus the EP populations leading to the strongest resonance can be identified with respect to the EP energy. In the following, the discussion of the AE stability focuses on the modes with the largest growth rates, that is to say, the dominant AEs.

Simulation results
This section is dedicated to analyzing the AE stability in different QPS configurations and NBI operational regimes.

Two field period configuration and vacuum equilibria
First, the AE stability is analyzed in the two period configuration and vacuum equilibria with respect to the NBI operational regime. Figure 2 shows the normalized growth rate (γτ A0 ) and frequency of the helical families n = 1, 3 and n = 2, 4 dominant modes as a function of the EP β and energy. Table 2 compares simulations with the same EP density although different EP energies. The maximum growth rates and the strongest resonance is observed if T f = 90 keV for n = 1, 3. The AEs with the largest growth rates show a frequency above 120 kHz. Figure 3 shows the AE growth rate and frequency with respect to the EP β (fixed T f = 90 keV) and T f (fixed EP . Graph (a) shows an increase of the AE growth rate as the EP density increases (equivalent to a stronger NBI injection power). The critical EP β to destabilize the AEs is 0.001. Panels (b) and (d) indicate the growth rate and frequency of the AEs destabilized by EP populations with different energy. The AE growth rate decreases as the EP energy increases at fixed EP β because the EP density decreases proportionally to the EP energy. The maximum growth rate is observed if the EP energy is 55 keV for the n = 1, 3 helical family and 30 keV for the n = 2, 4 helical family.
Next, each toroidal family is studied separately. Figure 4 shows the growth rate and frequency of the AEs triggered by the toroidal families n = 1 to 4 for different values of the EP β and energy. The AEs triggered by the n = 1 and n = 2 toroidal families show almost identical growth rate, frequency and eigenfunction structure compared to the AEs triggered by the n = 1, 3 and n = 2, 4 helical families (figure 2), respectively. Consequently, the helical couplings have a minor impact in the simulations. Table 3 shows the growth rate of the n = 1 to 4 toroidal families with respect to the EP energy (fixed the EP density). The AEs with the largest growth rate are triggered by the n = 1 and 2 toroidal families.  Figure 5 indicates the AE growth rate and frequency with respect to the EP β (fixed T f = 90 keV) and T f (fixed EP β = 0.03) for n = 1 to 4 toroidal families. The frequency jumps observed as the EP β increases indicate a change of the dominant mode between different AE families (panel c), for example, there is a transition between a dominant EAE to a NAE for the n = 4 toroidal family if the EP β is larger than 0.07. The AE stability trends of the n = 1 and n = 2 toroidal families are almost the same with respect to the n = 1, 3 and n = 2, 4 helical families profiles (figure 3), respectively. Nevertheless, the growth rate of the AEs triggered by the n = 1 and 2 toroidal families is 10%-20% smaller compared to the AEs destabilized by the n = 1, 3 and n = 2, 4 helical families. Consequently, the effect of the helical couplings in the AE stability is rather small. Figure 6, panel (a), shows the Alfvén gaps of the toroidal families 1 to 6 as well as the radial location and frequency range of the AEs triggered by the n = 1, 3 helical family. The TAE gap is small due to the symmetry of the magnetic field (approximately uniform in the poloidal direction). However, the residual poloidal asymmetry in the magnetic field also induces a rather wide EAE gap because shear Alfvén waves with m and m + 2 poloidal mode numbers cause a destructive interaction for a wider range of frequencies. This is caused by the fact that in QPS the poloidal symmetry is not exact and deviations from QP symmetry are present that increase towards the plasma edge. Specifically, the QP-optimization method gave a stronger priority to reducing the m = 1, n = 0 component of Bmn than to the higher m's (for example, the m = 2, n = 0 Bmn). As a result, significant m = 2 components remain in the magnetic field. This is an example as to how 3D plasma shaping can influence Alfvén gap structures and motivates future research on this topic (the QPS optimization did not include a target related to Alfvén gap widths). It should be recalled that the β-induced Alfvén eigenmodes (BAE) gap appears in panel (a) because the effect of the sound wave is included in the continuum calculation even though the analysis is dedicated to a vacuum case. The BAE gap shows an upper frequency bound of 25 kHz. There is a slender TAE gap that becomes wider towards the periphery of the plasma (r/a > 0.65) around 60-120 kHz. Also, there is a wide EAE gap all along the plasma radius, between 60-160 kHz in the inner-middle plasma and above 125 kHz at the plasma periphery. Small NAE gaps appear above 200 kHz. The

Three field period configuration and vacuum equilibria
The periodicity of QPS magnetic trap may affect the AE stability. This section is dedicated to the analysis of the AE stability if the QPS magnetic field has three periods.
The equilibrium modes in the three period configurations are n = 0, N, 2 * N, . . . with N the device magnetic field period. The helical couplings included in the simulation for a n = ±1 perturbation are n = 0 ± 1 = 1 and −1 as well as n = 3 ± 1 = 2 and 4. Thus, the helical mode family is n = 1 − 2 − 4. Figure 7 shows the growth rate and frequency of the AEs triggered by the n = 1, 2, 4 helical family for different values of the EP β and energy. The averaged growth rate is 80% smaller compared to the 2P Vacuum case, pointing out the AE stability improves if the number of field period increases to three. Table 4 shows the AE growth rate for simulations performed with the same EP density and different T f values for the n = 1, 2, 4 helical family. The AE with the largest   growth rate is induced by the EP population with T f = 180 keV showing a frequency of 240 kHz. If EP population energy is T f = 90 keV, AE with a frequency of 169 kHz is triggered. Figure 8 shows the AE stability trends regarding the EP β (fixed T f = 90 keV) and T f (fixed EP β = 0.03) triggered by the helical family n = 1, 2, 4. Panels (a) and (c) indicate AEs  Figure 9 shows the growth rate and frequency of the AEs triggered by the n = 1 to 4 toroidal mode families for different values of the EP β and energy. The AE stability trends obtained for the n = 2 toroidal family are similar to the trends obtained for the n = 1, 2, 4 helical family. Table 5 shows the AE growth rate in simulations with the same EP density and different EP energy. The largest AE growth rate is obtained for the EP population with T f = 180 keV triggering an EAE with a frequency around the 240 kHz destabilized between the middle-outer plasma region. No dominant HAEs are observed in the simulations. Figure 10 shows the AE growth rate and frequency with respect to the EP β (fixed T f = 90 keV) and T f (fixed EP β = 0.03) for n = 1 to 4 toroidal families. Panels (a) and (c)    Summarizing, the simulations for the 3P Vacuum case show AEs with a growth rate 1/5 compared to the 2P Vacuum case and an EP β threshold of 0.001 for all cases except for n = 1 toroidal family, with threshold of β = 0.05. This is explained by the slender continuum gap in the middle-outer plasma region due to a mixed TAE/HAE gap and an EAE gap displaced at a high frequency range in the middle-outer plasma  region. The unstable AEs observed in the simulations are TAEs and EAEs triggered in the middle plasma region with frequencies between 100-185 kHz. Despite the presence of HAE gaps in the plasma periphery, no dominant HAEs are identified in the simulations.

Two field period configuration and finite β equilibria
This section is dedicated to calculating the AE stability in QPS configuration with two periods and finite β, providing an improved description of the Alfvénic activity during the device operation. The study includes the analysis of the sub-dominant modes, NBI deposition region and FLR stabilization effects. Figure 12 shows the growth rate and frequency of the AEs triggered by the helical families n = 1, 3, 5 and 2, 4, 6 for different EP β and energy values. The EP β threshold is 0.06 for the n = 1, 3, 5 helical family and 0.03 for the n = 2, 4, 6 helical family, much higher compared to the vacuum cases. The AEs frequency triggered above the EP β threshold is similar to the vacuum cases, although the growth rate is lower. Therefore, finite β case is more stable with respect to the vacuum case. Table 6 shows the AE growth rate in simulations with the same EP density and different EP energies. The strongest resonance is observed for the EP population with T f = 90 keV and the n = 1, 3, 5 helical family triggers AEs with larger growth rate compared to the n = 2, 4, 6 helical family. Figure 13 shows the AE stability trends regarding the EP β (fixed T f = 90 keV) and T f (fixed EP β = 0.07) triggered by the n = 1, 3, 5 and n = 2, 4, 6 helical families. Panel (a) and (c) show the EP β threshold is 0.07 for the n = 1, 3, 5 helical family and 0.05 for the 2, 4, 6 helical family. Panels (b) and (d) indicate that the AEs are stable if T f > 130 keV for both helical families. Figure 14 shows the growth rate and frequency of the AEs triggered by the n = 1 to 5 toroidal families for different EP β and energy values. The AEs growth rate is lower with respect to the vacuum case and the n > 1 toroidal mode families mainly trigger NAEs with frequencies above 400 kHz. Table 7 shows the AE growth rate in simulations with the same EP density and different EP energies. The strongest resonance is caused by the EP population with T f = 130 keV for the n = 1, 2 and 3 toroidal mode families and by the EP population with T f = 90 keV for the n = 4 and 5 toroidal mode families. Figure 15 shows the AE stability trends regarding the EP β (fixed T f = 90 keV) and T f (fixed EP β = 0.07) triggered by the n = 1 to 5 toroidal families. Panels (a) and (c) indicate the β EP threshold to destabilize n = 1 AEs is 0.001, 0.02 for the 2 and 3 AEs, 0.03 for the n = 4 AEs and 0.06 for the n = 5 AEs. Panels (b) and (d) show that n = 3 AEs are only triggered by EP populations with energies below 130 keV as well as n = 4 and 5 AEs by EP populations with an energy below 90 keV.    transitional mode between a 1/3 − 1/4 TAE and a 1/3 − 3/9 HAE with 98 kHz triggered by the EP population with 55 keV. Panel (g) shows an unstable 2/6 − 2/7 TAE with 113 kHz. All the AEs are triggered in the middle-outer plasma region.
The growth rate of the sub-dominant modes (modes with a lower growth rate than the dominant mode) is calculated for the n = 1 toroidal family if T f = 90 keV and β = 0.02. Table 8 shows the growth rate of the sub-dominant modes identified. Only sub-dominant n = 1 AEs are identified, showing growth rates one and two orders of magnitude lower with respect to that of the dominant mode. Thus, the effect of the sub-dominant modes on the plasma stability should be negligible. Figure 17 shows the eigenfunctions of the sub-dominant modes, 1/3 − 1/4 TAEs with f = 80-69 kHz destabilized between the middle-outer plasma region. Now, the effect of the NBI deposition region on the AE stability is analyzed for n = 1 toroidal family if T f = 90 keV and β = 0.02. This study is performed displacing the gradient of the EP density profile from the magnetic axis (on-axis case) to the plasma periphery (off-axis cases). The profiles used in the analysis are indicated in the figure 18, panel (a). The parameter r peak indicates the radial location of the EP density profile gradient. Panels (b) and (c) indicate the growth rate and frequency of n = 1 to 5 dominant AEs if the NBI radial deposition region is modified. Simulations for an on axis NBI injection only leads to the destabilization of n = 1 AEs, although simulations with off axis NBI injection leads to the destabilization of n = 2 and 3 AEs if r peak = 0.5 as well as n = 4 and 5 AEs if r peak = 0.7. Thus, an off axis NBI deposition causes a stronger destabilization of the AEs with respect to the on axis NBI injection. This is explained by the wide TAE gap located in the plasma periphery. Figure 19,  It should be noted that the EP β decreases as the NBI is deposited off-axis because the slowing down time of the EP is shorter compared to on-axis configurations. The effect of the EP β decrease as the NBI deposition is displaced off-axis is not included in the study because there is not a straight forward relation of the EP β and the NBI deposition region, thus the EP β is assumed constant for simplicity. Consequently, a further analysis is required to identify the optimal configuration with respect to the NBI injection intensity and radial deposition region to improve the AE stability. Nevertheless, the present study already provides useful information because the AE stability strongly depends on the radial location of EP density profile gradient with respect to the continuum gaps.
In summary, the two field period finite β case shows an improved AE stability compared to the vacuum cases because the main part of the unstable modes are high frequency modes strongly damped by the EP FLR effect. The AEs triggered by the helical family 1, 3, 5 show a growth rate 10%-50% smaller compared to the vacuum case. Regarding the AEs triggered by the helical family 2, 4, 6, the growth rate is 20% if the AE is triggered by EP populations with energies below 90 keV. On the other hand, the AEs triggered by EP populations with an energy of 90 keV or higher, the AE growth rate is similar or slightly higher compared to the vacuum case simulations. The growth rate of the sub-dominant modes is one order of magnitude smaller with respect to the dominant modes and only n = 1 sub-dominant modes are identified. Off-axis NBI injection leads to a stronger destabilization of the AEs with respect to the on-axis NBI injection because there is a wide TAE gap in the plasma periphery.

Conclusions
The gyro-fluid FAR3d code is applied to calculate the AE stability in different QPS configurations. This study belongs to a compilation of studies dedicated to analyzing the heating efficiency in different magnetic field configurations, including tokamaks as well as stellarators that takes advantage of different quasi-symmetries.
The simulations indicate the destabilization of TAE/EAEs by the 1, 3 and 2, 4 helical families in the frequency range of 90-370 kHz (AEs with the largest growth rates) for the configuration with two magnetic field periods for vacuum equilibria. The EP β threshold is 0.001 and the strongest resonance is caused by the EP population with 55 keV.
The AE stability is improved in the configuration with three magnetic field periods, showing a growth rate five times smaller than that in the two field period configuration. The improvement of the AE stability is linked to the presence of slender continuum gaps between the middle-outer plasma region, showing a mixed TAE/HAE gap. More simulations are required to identify the optimal QPS configuration with respect to the number of magnetic field periods, although an improvement in AE stability is observed if the number of field periods is increased from two to three.
The simulations for the finite β equilibria in the two field period case indicate the destabilization of TAE/EAEs by the helical families n = 1, 3, 5 and n = 2, 4, 6. The EP β threshold is higher compared to the vacuum cases, above the EP β = 0.03. The destabilization of a 1/3 − 3/9 HAE with f = 98 kHz is identified in the simulations at the plasma periphery triggered by the EP population with an energy of 55 keV.
The analysis indicates on-axis NBI injection improves the AE stability with respect to off-axis NBI injection. If the NBI is deposited at r/a ⩾ 0.5, AEs are destabilized in the middleouter plasma region. The improved stability of the on-axis NBI injection is linked to a slender TAE gaps in the inner plasma as well as rather high EP β threshold of the EAEs.
QPS configurations with two periods and finite β as well as the vacuum case with the three periods show the presence of helical gaps in the middle-outer region of the plasma, but not in the vacuum two period configuration. Nevertheless, the simulations indicate TAEs and EAEs show the largest growth rate and only dominant HAEs are identified in the configuration with two periods and finite β.
The analysis indicates the QPS configurations with three periods for an equilibria with a finite β may show a robust AE stability and an efficient heating by a hypothetical NBI with the same properties as the tangential beams of the LHD device, particularly if the NBI is deposited on-axis. If the NBI power injection leads to an EP β below 0.03, the AEs are stable. If the EP β threshold is exceeded, n = 1 and 2 TAE/EAE can be triggered in the middle plasma region.