MHD stability trends and improved performance of LHD inward-shifted configurations: The role of the neutral beam current drive and thermal plasma density Open Access

Several actuators can be applied to modify the MHD stability of nuclear fusion plasmas. For example, the neutral beam (NBI) in nuclear fusion devices can alter the magnetic field topology by the generation of non-inductive currents.^1,2 Neutral beam current drive (NBCD) can be used to achieve steady-state operation in advanced tokamaks^3–5 or to modify the magnetic field configuration.^6–9 NBCD can also affect the stability of pressure gradient-driven modes (PGDM)^10–15 and Alfvén eigenmodes (AE).^16,17

NBCD generation was analyzed in the Large Helical Device (LHD) experiments,^18,19 leading to the stabilization of toroidal and global Alfvén eigenmodes (TAE/GAE),²⁰ as well as PGDM.^21,22 Also, theoretical analysis identified scenarios with an improved MHD stability with respect to the NBCD intensity.²³ The destabilization of PGDM deteriorates LHD plasma confinement^24–28 and unstable AEs reduce the plasma heating efficiency,^20,29–36 leading to a decrease in the device performance.^37–46 

Another actuator affecting the MHD stability is the thermal plasma density and temperature. PGDM stability is closely linked to the thermal plasma β of the discharge, triggered once a threshold of the pressure gradient is exceed.²⁴ Likewise, energetic particle modes (EPM) as the energetic-ion-driven resistive interchange mode (EIC) are stabilized by off-axis ECH injection, increasing the temperature of the plasma periphery.^47–49 

The operational regime of the NBI is important to improve the AE stability too. Optimized DIII-D discharges with lower AE activity were observed balancing the NBI voltage and power injection while keeping constant the EP β⁵⁰ or by displacing outward the NBI deposition region.^51,52

LHD plasma is heated by three tangential NBI lines parallel to the magnetic axis injected at the plasma core with a beam energy of 180 keV as well as two NBIs perpendicular to the magnetic axis injected in the plasma periphery with a beam energy of 32 keV.^53–55 An unbalanced injection of the tangential NBIs leads to the generation of a net plasma current, modifying the rotational transform profile as well as the stability of PGDM and AE.

The present study is dedicated to identify operation scenarios of LHD inward-shifted configuration showing an improved MHD stability with respect to the NBCD, thermal plasma density, and NBI operational regime. On that aim, several experiments in the 22th and 23th LHD campaigns were performed exploring the effect of the NBCD and thermal plasma density on the PGDM and the AE stability, testing unbalanced NBI heating patterns leading to a net plasma current up to 60 kA/T. In addition, the linear stability of AEs is analyzed using the gyro-fluid code FAR3d to reproduce the stability trends observed in the experiments.^56–59 This study explores for the first time possible operation scenarios showing, at the same time, an improved thermal plasma confinement due to a reduction of the PGDM activity as well as an enhanced plasma heating by minimizing the AE activity and mandatory requirements for the economical viability of future fusion reactors.

The FAR3d code solves the reduced nonlinear resistive MHD equations coupled with the moment equations of the energetic ion density and parallel velocity.^60–62 The Landau damping/growth is included by Landau closure relations, reproducing the linear wave–particle resonance effects on six field variables that evolve from three-dimensional equilibria generated by the VMEC code.⁶³ 

This paper is organized as follows: A description of the experiments performed is done in Sec. II. The numerical scheme and equilibrium properties are described in Sec. III. The analysis of the PGDM and the AE stability with respect to the NBCD is done in Sec. IV. The AE stability for different thermal plasma densities is studied in Sec. V. The effect of the NBI operational regime on the AE stability is discussed in Sec. VI. Improved PGDM and AE stability in the discharge 167 800 by the control of NBCD thermal plasma density is discussed in Sec. VII. Next, the analysis results are discussed in Sec. VIII. Finally, the conclusions of this paper are presented in Sec. IX.

Experiments in the 22th and 23th LHD campaign, shot numbers 167 777–167 808 and 176 665 to 176 676, explored LHD inward-shifted configurations (counter-clock-wise orientation of the magnetic field) with co-, counter-, and balanced NBCD for different thermal plasma densities, identifying the operation scenarios leading to the destabilization of PGDM and AE in a hydrogen plasma.

The analysis of the MHD stability for the discharge 167 777, shown in Fig. 1, indicates the destabilization of $n/m=1/2$ PGDM [panels (c) and (d)] along the ctr-NBCD phase (NBI 2 injection) above a given threshold of the plasma toroidal current [I_p, panel (b)]. On the other hand, 1/2 PGDM is stable in the co-NBCD phase (NBI 1 and 3 injection) although 3/2, 4/3, and 1/1 PGDM are triggered at the plasma periphery. The averaged thermal β is rather constant during both NBCD phases [panel (a)], and thus, the evolution of the PGDM stability is caused by the NBCD, particularly due to the modification of the iota profile. NBI 1 and 3 induce a co-NBCD leading to an up-shift of the iota profile between the inner-middle plasma region. On the other hand, NBI 2 generates ctr-NBCD leading to a down-shift of the iota profile. Magnetic fluctuation diagnostics in LHD consist of 64 Mirnov coils and 24 saddle loops installed inside the vacuum vessel in order to measure magnetic fluctuations.

Figure 2 shows the AE activity in the discharge 167 794 along the ctr-NBCD phase and in the discharge 167 807 along the co-NBCD phase [panels (b) and (e), respectively]. In the discharge 167 794, if $Ip<−5$ kA/T, the AE around 60 kHz is stabilized and the AE around 45 kHz shows a frequency up-shifts to 55 kHz [panel (c)]. On the other hand, the discharge 167 807 shows the destabilization of AEs with a frequency around 50 kHz if I_p > 10 kA/T [panel (f)]. Again, the averaged thermal β (both thermal plasma density and temperature, data not shown) is rather constant, thus the modification of the AE stability is mainly caused by the effect of the NBCD.

In the following, a subgroup of experiments is selected to analyze the stability of PGDM and AE, identifying the mechanism leading to the mode stabilization/destabilization. The actuators in the MHD stability discussed in the study are the effect of the NBCD on the iota profile, the averaged thermal β (both thermal plasma density and temperature), and the NBI operational regime (NBI voltage, injection power, and radial deposition region). The evolution of the iota profile is measured using the motional stark effect (MSE) diagnostic, although MSE data are only available if NBI 3 is injected.

FAR3d code solves a reduced set of equations for high-aspect ratio configurations and moderate β-values (of the order of the inverse aspect ratio), retaining the toroidal angle dependency in an exact three-dimensional VMEC equilibrium.⁶⁴ The destabilizing effect of the EP particles is included toward moments of the gyro-kinetic equation, particularly the EP density (n_f) and the EP velocity parallel to the magnetic field lines (⁠ $v||f$⁠). The correct model calibration requires performing gyro-kinetic simulations to calculate the Landau closure coefficients in the gyro-fluids simulations, matching the analytic TAE growth rates of the two-pole approximation of the plasma dispersion function with a Lorentzian energy distribution function for the energetic particles. The lowest-order Lorentzian can be matched either to a Maxwellian or to a slowing-down distribution by choosing an equivalent average energy. Furthermore, details of the model equations can be found in Refs. 65 and 66. A Maxwellian equilibrium EP distribution is chosen, which has the same second moment, effective EP temperature, as the slowing-down distribution as defined in the  Appendix.

Not all the resonances identified by the simulations should lead to the destabilization of AEs in the experiment, because the drive is determined by the gradient of the phase space distribution and the gradient depends on the phase space shape of the distribution function of the EP. Nevertheless, this information is useful for future optimization studies. It is important to point out that the model reproduces the destabilizing effect of the passing EP; hence, the effect of highly anisotropic beams or ICRF-driven EP cannot be modeled by the present version of the code. However, the pitch angle of the EP generated by the tangential N-NBI in LHD plasma should be small. In addition, due to the Maxwellian distribution function used for the EP model, the co-EP (pitch angle is 0) and ctr-EP (pitch angle is π radians) lead to the same resonance (same frequency and growth rate) although the mode propagates in the opposite direction. Thus, the observed modes can be caused by co-passing EPs generated by the N-NBI.

Previous benchmarking studies validated FAR3d code results with respect to gyro-kinetic and hybrid codes.⁶⁷ FAR3d code was applied to model the MHD activity in LHD,^23,68–70 DIII-D,^50,71–74 TJ-II,^59,75,76 and Heliotron J^77–79 plasma, showing a reasonable agreement. Also, the nonlinear version of the code reproduced the sawtooth-like, internal collapse events, EIC burst, and MHD burst observed in LHD plasma^66,80–84 as well as the AE saturation in DIII-D plasma.⁸⁵ FAR3d code was applied to analyze the AE stability in plasma with multiple EP populations for LHD, DIII-D, ITER, and CFETR^86–88 and to predict the MHD stability of future devices as JT60SA and CFQS.^65,89

This section is dedicated to introduce the equilibrium used in the simulations. In addition, the simulation parameters are described.

A set of equilibria is calculated using the VMEC code⁶³ mimicking the iota profile down-shift in the ctr-NBCD phase of 167 805 discharge as well as the iota profile evolution in the discharge 167 800 at t = 5.0, t = 5.5, and 5.8 s. The electron density and temperature profiles are reconstructed using Thomson scattering and electron cyclotron emission data. The birth energy of the particles injected by the tangential NBI is $Tf(0)=180$ and 40 keV for the perpendicular NBI. The reference model assumes the nominal energy of the passing EP injected by the tangential NBI is $Tf(0)=100$ keV, constant along the normalized minor radius for simplicity. No plasma heating by the perpendicular NBIs is applied in the discharges analyzed.

Table I indicates the dynamic and equilibrium toroidal and poloidal modes included in the simulations. Equilibrium modes (n = 0) represent the magnetic trap configurations and do not evolve in time. Dynamic modes (⁠ $n=1,2,3,4$⁠) represent the perturbations evolving in time. The analysis of the PGDM stability includes n = 1 to 4 toroidal mode families. The mode selection adds all the resonant modes between the magnetic axis and the plasma periphery. The poloidal mode selection is extended in the simulations dedicated to study the AE stability as the iota profile is up/down-shifted, including all the rational surfaces between $ι=0.1$ and 1.5. The number of point of the radial grid is 1000 and the magnetic Lundquist number is $5×106$⁠. The eigensolver version of the code is used to calculate the stability of the dominant (largest growth rate) and sub-dominant modes. The analysis of the AE stability is limited to the n = 1 toroidal family because it triggers AEs with the largest growth rate. EP Finite Larmor radius (FLR) damping effect is included in all the simulations (EP Larmor radius is 0.023 m) expect for the analysis performed in Sec. VI.

This section is dedicated to analyze the effect of the NBCD on PGDM and AE stability. The deformation of the iota profile caused by the NBCD modifies the radial location of the rational surfaces and the structure of the continuum, thus PGDM and AE stability changes.

The deformation of the iota profile by the NBCD can lead to the transition of rational surfaces from resonant to non-resonant. In addition, the radial location of the resonant rational surfaces changes, modifying the PGDM stability. It should be noted that the EP effect on the stability of PGDM is not included in the analysis for simplicity.

Figure 3 shows 1/2 PGDM is stable along the co-NBCD phase of 167 800 discharge [panels (b)–(d)]. The stabilization of 1/2 PGDM at the beginning of the co-NBCD phase is caused by the distortion induced by the NBCD on the iota profile, because the maximum averaged thermal β is reached after the mode stabilization [panel (a)].

Figure 4 shows the evolution of the iota profile for the discharge times indicated in Fig. 3. The iota profile is measured using the MSE diagnostic at t = 5.5 and 5.8 s although the iota profile at 5.0 s is calculated by the VMEC code (no MSE data available) using a simplified current profile $Ip=16.8(1−s)$ kA (consistent with the Rogowski coil data). The iota profile increases along the co-NBCD phase after the screening current induced by the NBIs 1 and 3 begins to decay [panel (a)]. Consequently, the up-shift of the iota profile changes the radial location of the rational surfaces, leading to the rational surfaces 1/2 to be non-resonant at t = 5.8 s and explaining the mode stabilization.

In summary, co-NBCD leads to an improvement of the device performance by stabilizing 1/2 PGDM, inducing a transition of the 1/2 rational surface from resonant to non-resonant. This result is consistent with previous theoretical studies dedicated to analyze the effect of the co-NBCD intensity on the PGDM stability.²³ It must be recalled that the rather large uncertainty of MSE measurements requires a validation of the analysis conclusion once more accurate iota profile measurement are available. Nevertheless, the PGDM stability trends with respect to the NBCD evolution along the experiment are rather robust.

The modification of the iota profile by the NBCD leads to a reconfiguration of the continuum gap structure. Consequently, if the width, radial location, and frequency range of the gaps change, the effect of the continuum damping on the AE stability is also modified leading to the stabilization or further destabilization of the AEs.

Figure 5 shows the stabilization of the AE activity observed during NBI 2 injection once the ctr-NBCD intensity is above 14 kA/T from t = 4.8 s [panels (c)]. PGDM are stable because the averaged thermal $β<0.3%$ is below the PGDM destabilization threshold [panel (a)]. The analysis is performed during the ctr-NBCD phase so no MSE data are available to track the evolution of the iota profile [panels (b) and (d)]. Nevertheless, the MSE data at the beginning of the co-NBCD phase indicate the iota profile may reach values below 0.3 during the ctr-NBCD phase.

The location of the gradient is controlled by the parameter r_peak (equal to 0.25) and the flatness by δ_r (equal to 3). It should be noted that no calculations of neutral beam deposition are performed for this discharge; instead, a set of parameter studies are done with respect to the EP density profile and β to approximate the resonance observed in the experiment.

Figure 7(a) shows the continuum gaps calculated for different iota profiles using the code Stellgap including the effect of the sound wave.⁹³ If extreme cases are compared, the down-shifted iota profile with $ιmin=0.1$ (black line) and the up-shifted with $ιmin=0.5$ (pink line), there are important differences between the continuum gaps. The $ιmin=0.1$ case shows slender gaps in frequency range and radial width, mainly located between the middle-outer plasma regions. Also, the frequency range of TAE and EAE gaps in the middle plasma is rather low, between $70−90$ and $125−140$ kHz, respectively. Thus, AEs triggered in this configuration may have narrow eigenfunctions and experience an enhanced stabilizing effect by continuum and FLR dampings (Alfvén gaps are partially closed and there is a radial up-shift of the continuum frequency range). Consequently, if the NBI injection is performed on-axis and the EP density gradient is located in the inner plasma region, this configuration is particularly robust with respect to the AE stability. On the other hand, $ιmin=0.5$ case shows wide gaps in frequency range and radial width. The TAE gap is displaced to a high-frequency range, between $200 and 250$ kHz, thus EP FLR effect may have a large stabilizing effect on the TAEs although the continuum damping is rather weak in this configuration, particularly in the inner plasma region. Also, unstable TAEs may show wide eigenfunctions that can spread along all the minor radii. Regarding intermediate configurations, $ιmin=0.3$ case (green line) shows a TAE gap in the inner plasma region between the frequencies $60−140$ kHz, as well as a wide BAE gap near the magnetic axis. Figure 7(b) indicates the growth rate and frequency of the dominant n = 1 AEs in the configuration $ιmin=0.3$ for different EP β values. The EP β threshold for AE destabilization is 0.002. Growth rate and frequency are larger as the EP β increases. Figure 7(c) indicates the growth rate and frequency of the dominant n = 1 AE in configurations with different iota profiles (EP $β=0.01$⁠). Only n = 1 family is shown because the growth rate is $50%$ higher compared to the dominant n = 2 AE. The simulations indicate the configurations with the lowest AE growth rate are the $ιmin=0.4$ and 0.1 cases. The decrease in the AE growth rate in $ιmin=0.1$ configuration is linked to the narrow continuum gaps, although for the $ιmin=0.4$ case, it is caused by the continuum band that appears the inner plasma where the EP density profile gradient is located.

Figure 8(a) compares experimental data with simulation results. The magnetic probe (MP) data indicate a decrease in the AE frequency from 65 to 60 kHz between t = 3.5 and 4.0 s although it increases from 60 to 70 kHz between t = 4.0 and 4.7 s before the AE stabilization. Panels (b) and (c) show the measured toroidal and poloidal mode numbers of the instability, respectively. The instability is an n = 1 AE with a dominant m = 1–2. The vacuum equilibria calculated by VMEC code shows a $ιmin≈0.4$⁠; thus, the iota profile up-shift caused by screening currents at the beginning of the ctr-NBCD phase may lead to a $ιmin>0.4$⁠. The up-shift of the iota profile is approximately reproduced by $ιmin=0.5$ configuration, showing a dominant 1/2 Global AE (GAE) with f = 66 kHz triggered nearby the local minima of the continuum gap, close to the magnetic axis [panel (d)]. It should be noted that this AE is triggered inside the TAE gap, showing a weak coupling between the 1/2 and 1/3 modes. That means, this AE is a transitional mode between a 1/2 GAE and a $1/2−1/3$ TAE. Transitional modes are observed as the simulation result if the model is unable to distinguish between two different modes located at a similar frequency range and radial location. Once the effect of the screening current vanishes, the iota profile begins the down-shift along the ctr-NBCD phase. The frequency of the dominant AE calculated in the simulation as the iota profile down-shift may agree with the AE frequency sweeping observed in the MP data. The simulation for $ιmin=0.4$ case indicates a decrease in the AE frequency from 65 to 50 kHz that corresponds to a 1/2 β induced AE (BAE) triggered in the inner plasma region [panel (e)]. The frequency jump obtained in the simulations is 10 kHz larger compared to the experiment, predicting the destabilization of AEs at different gaps. On the other hand, the measured instability frequency shows a continuous variation, that is to say, there is no transition between AEs at different gaps. This inconsistency that can be explained by the uncertainty in the thermal plasma and EP profiles used in the study. Also, the EP density profile evolves along the experiment due to the outward EP transport induced by the AEs, leading to a decrease in the profile gradient although an increase in the profile flatness. This effect is not captured in the simulations because the EP density profile is fixed. Hence, the simulations show the same AE frequency trends with respect to the experiment as the iota profile evolves, although the dominant mode identified is incorrect. Analyzing the continuum structure, the instability could be a transitional 1/2 GAE to $1/2−1/3$ TAE located in the middle plasma, nearby the minimum of the continuum at $r/a=0.5$⁠. The unstable modes in the configurations with $ιmin=0.5$ and 0.4 do not intersect the continuum at $r/a<0.1$⁠, panels (h) and (i), reason why these modes are identified as gap modes (the continuum data are not shown at $r/a<0.1$⁠, the radial location of the first flux surface of the equilibrium, although the smooth radial variation of the continuum frequency range allows the extrapolation). The simulations for $ιmin=0.3$ and 0.2 configurations indicate an increase in the AE frequency to 65 and 72 kHz, respectively. The dominant AE is a $1/2−1/3$ TAE in $ιmin=0.3$ configuration [panel (f)] and a $1/3−1/4$ TAE in 0.2 case [panel (g)], both triggered in the inner plasma region. The eigenfunction of the AE in $ιmin=0.1$ configuration is not included in the figure, and it is a $1/5−1/7$ EAE with f = 89 kHz triggered in the inner plasma region. The AE identification is done considering the radial location and frequency range of the AEs with respect to the continuum [panels (h)–(k)]. The MP data show a faint AE activity from 4.7 s that reaches 80 kHz at t = 5.0 s, indicating the AE stabilization by the enhancement of the continuum and EP FLR damping effects as the iota profile down-shifts. It should be noted that the simulations do not reproduce the measured AE toroidal mode number, and the predicted poloidal mode number is smaller. Such discrepancy is caused by an underestimation of the destabilizing effect of the m = 1 mode in the simulations, probably linked to the inaccurate iota profile used in the simulations compared to the experiment.

In summary, the effect of the NBCD on the AE stability is linked to the modification of the continuum gaps as the iota profile changes. The ctr-NBCD leads to configurations with slender gaps in frequency range and radial width, particularly at the inner plasma region; thus, continuum and FLR damping effects are stronger and the EP β threshold for the AE destabilization is higher. Nevertheless, TAE and EAE gaps at the middle-outer plasma that can lead to the destabilization of AEs if there is a EP density profile gradient in this plasma region. Co-NBCD configurations show wide continuum gaps at high-frequency ranges; for example, the TAE gap is observed around 200 kHz covering all the minor radius. Consequently, co-NBCD configurations have a stronger EP FLR damping although weaker continuum damping compared to ctr-NBCD cases. Co-NBCD configurations also show a higher EP β threshold for AE destabilization although, if the drive is large enough to trigger an AE, the instability can cover main part of the minor radius. In addition, GAEs/BAEs can be destabilized nearby the lower bound of the TAE gap frequency range in the inner plasma region. In short, configurations showing an improved AE stability are identified above a given threshold of co- and ctr-NBCD intensity. These results are consistent with previous theoretical studies dedicated to analyze the effect of the co-NBCD intensity on the AE stability.²³ 

This section is dedicated to study the effect of the thermal plasma density on PGDM and AE stability. An increment of the thermal plasma density leads to a higher thermal β, eventually exceeding the PGDM stability limit as the pressure gradient increases.^24,28,94 Above the PGDM stability limit inward-shifted LHD configurations show unstable interchange modes.^25,27,41,95 A modification of the thermal ion density also changes the plasma Alfvén velocity, defined as $VA=B/Nμ0mpni$ with B the magnetic field intensity, μ₀ the vacuum magnetic permeability, n_i the thermal ion density, m_p the proton mass, and N the ion species mass/proton ratio. The resonance induced by the EP can be defined with respect to the ratio between the velocity of the EP $Vf=Tf/Nmpni$ and V_A. If $Vf/VA≈1$⁠, the EP resonance is strong and the EP β threshold to destabilize the AEs is lower. In addition, the EP β also changes as the thermal plasma density increases because the NBI deposition rate and the EP slowing down depend on the thermal plasma density too. Particularly, the NBI deposition rate (⁠ $PNBI,dep$⁠) increases linearly with the electron plasma density (⁠ $⟨ne⟩$⁠) for low-density plasma with $⟨ne⟩<1019$ m⁻³ saturating for higher plasma densities. On the other hand, the EP slowing-down time (τ_ep) is proportional to $Te3/2/ne$ (T_e is the electron temperature). Thus, the EP β dependency with the thermal plasma density is, roughly speaking, proportional to $PNBI,dep*τep$⁠. That means, configurations with an improve stability of PGDM and AE with respect to the thermal plasma density may not be the same. A large thermal plasma density enhances the PGDM activity although the decrease in V_A may lead to a weaker or a stronger EP resonance depending on the $Vf/VA$ ratio. The simulations performed include EP distribution functions modified with respect to the variation of the $Vf/VA$ ratio and EP β as the thermal plasma density change, re-scaling the EP density, and energy profiles to provide an approximated description of the AE stability trends with respect to the thermal plasma density.

Figure 9 shows the MHD activity in a set of discharges as the averaged thermal electron density (⁠ $⟨ne⟩$⁠) increases from 0.2 to $4.1× 1019$ m⁻³, leading to an enhancement of the averaged thermal β (β_th) from 0.21% to $1.65%$ and an increase in $Vf/VA$ ratio from 0.14 to 0.66 (sub-Alfvénic beam). Low-density discharges show AE activity in frequency ranges above 40 kHz (the yellow dashed horizontal line indicates the upper frequency limit of the BAE gap). The MHD activity observed in the discharges with $⟨ne⟩=0.2$ to $0.8×1019$ m⁻³ is Alfvénic [panels (a)–(c)], although AE and PGDM coexist if $⟨ne⟩=1.2×1019$ m⁻³ (⁠ $βth=0.65$⁠) once the PGDM stability limit is exceeded [panel (d)]. Keep increasing $⟨ne⟩$ causes the stabilization of the AEs although a further enhancement of the PGDM. Weak AE signal is still observed for $⟨ne⟩=1.7$ to $2.1×1019$ m⁻³ as the PGDM signal increases below 10 kHz due to the destabilization of 3/2, 1/1, and 4/3 interchange modes at the plasma periphery [panels (e) and (f)]. If $⟨ne⟩≥2.7×1019$ m⁻³, the AEs are stable and the PGDM activity is rather large [panels (g)–(j)], particularly in the discharge 167 785 that shows unstable 3/2, 4/3, and 1/2 interchange modes inducing large amplitude magnetic perturbations. Consequently, the experimental observations indicate PGDM and AEs optimization trends with respect to $⟨ne⟩$ are opposite. In addition, there are operation scenarios handicapped by both instabilities as the same time that should be avoided. High-density operations may show an improved plasma heating efficiency with stable AEs, although strong PGDM activity may deteriorate the plasma confinement.

Figure 10 shows the variation of the continuum gaps as the thermal plasma density increases from 0.2 to $2×1019$ m⁻³, indicating slender gaps in configurations with larger thermal plasma density. Figures 11(a) and 11(b) indicate the growth rate of the dominant n = 1 AE or EPM (black dots) and 1/2 PGDM (red dots) with respect to the thermal plasma density. In this simulation, $Vf/VA$ ratio changes because V_A is different for each thermal plasma density tested although V_f is fixed. The n = 1 AE with f  > 45 kHz is stable if $n≥1.6×1019$ m⁻³ consistent with the experiment observations. The sub-dominant mode with the largest growth rate is the n = 1 AE if $n<0.8×1019$ m⁻³, replaced by a 1/2 PGDM with $f≈25$ kHz for $n≥0.8×1019$ m⁻³. The scaling of the n = 1 AE/EPM growth rate is nonlinear with respect to the thermal plasma density, $γτA0=0.10−0.03n2.5$⁠, linked to the nonlinear dependency of the EP resonance intensity with the thermal plasma density, proportional to the $n$⁠. It should be noted that the enhancement of the EP destabilizing effect linked to the increase in the $Vf/VA$ ratio from 0.14 to 0.65 is compensated by the combination of a lower EP drive (EP β decreases as the thermal plasma density increases if $⟨ne⟩<1019$ m⁻³) and more radially localized AE/EPM eigenfunction as the Alfvén gaps gets slender (stronger FLR damping effect), leading to an improved AE/EPM stability. The increment of the n = 1 AE/EPM frequency obtained in the simulations as the thermal plasma density increases is counter-intuitive because the AE frequency should decrease as the thermal plasma increases because $fA0=1/τA0=B0/(R2μ0min)1/2$⁠. Nevertheless, it must be clarified that the growth rate and frequency indicated in the figures correspond to the dominant AE or EPM identified in the simulation. That means, the combined variation of the Alfvén continuum structure, EP resonance properties, EP β, and FLR effects as the thermal density increases, leads to the destabilization of AE or EPM with different characteristics. In particular, the large increase in the AE frequency if $n=1.2–1.4×1019$ m⁻³ is caused by the destabilization of a transitional mode between a 1/3 EPM and an AE in the EAE and NAE gaps. The scaling of the 1/2 PGDM growth rate with the thermal plasma density is almost lineal, $γτA0=0.08+0.03n1.07$⁠. The green stars in panels (a) and (b) show the 1/2 PGDM growth rate and frequency if the decrease in the EP β as the thermal ion density increases is included in the simulations. Decrease the EP β in the simulation is equivalent to a proportional down scale the EP density profile [Fig. 6(b)], that is to say, reducing the EP population in the plasma. The growth rate slightly increases because the EP population induces a stabilizing effect on the 1/2 PGDM, smaller as the EP β decreases for a larger thermal ion density. In addition, 1/2 PGDM frequency decreases from f = 25 to 6 kHz, indicating the weaker coupling between the shear Aflvenic wave with the PGDM as the EP β decreases. Figures 11(c) and 11(d) show the eigenfunction of the n = 1 AE and 1/2 PGDM for $n=0.6×1019$ m⁻³. The n = 1 AE is a transitional mode between a 1/3 GAE and a $1/2−1/3$ TAE destabilized in the inner plasma region and the 1/2 PGDM is triggered in the middle plasma region. The 1/3 GAE or $1/2−1/3$ TAE correspond to the AE observed in the discharges with $n≤2.1×1019$ m^{− 3}. The 1/2 PGDM is consistent with the MHD activity measured in the frequency range between 20 and 40 kHz in the discharges with n = 1.7 to $2.7×1019$ m⁻³. It must be mentioned that no PGDM activity with f > 10 kHz is measured in the discharges with $n>2.7×1019$ m⁻³, consistent with the 1/2 PGDM frequency down sweeping predicted in the simulations.

In summary, the stability of PGDM and AEs shows opposite trends with respect to the thermal plasma density. AEs are stable at high thermal plasma densities (above $2×1019$ m⁻³), because the continuum gaps are slender in radial width and frequency range, displaced at a lower frequency range. On the other hand, the increase in the thermal β as the thermal plasma density increases leads to the destabilization of PGDM as soon as the plasma pressure exceeds the PGDM stability limit. In addition, the MHD activity measured between 20 and 40 kHz is identified as an 1/2 PGDM. These modes are only observed if the thermal plasma density is in the range of 1.7– $2.7×1019$ m⁻³. Operations with a large thermal plasma density do not show MHD activity above 10 kHz, explained by the PGDM frequency down sweeping observed in the simulations as the EP β decreases due to the thermal ion density increase and the coupling of the PGDM with the shear Aflvenic wave is weakened.

This section is dedicated to analyze the effect of the NBI operational regime on the AE stability. The study includes the AE stability trends with respect to the NBI injection power that is linked to the amount of EP in the plasma (EP β), the NBI voltage associated with the energy of the EP (T_f), and the NBI deposition region (on-axis vs off-axis injection). In the simulations, modifying the EP β corresponds to an up/down scaling of the EP density profile, a change of the EP energy to an up/down scaling of the EP temperature profile (assumed constant along the minor radius for simplicity) and the NBI deposition region is approximated by a displacement of the EP density gradient radial location. FLR effects are not included in the simulations for simplicity.

Figure 12(a) shows the destabilization of AEs in the frequency range of 75 kHz during the ctr-NBCD phase of the discharge 167 800, after the 1/2 PGMD is stabilized by the iota profile up-shift. The analysis of theoretical NBI operational regimes is performed at the discharge time t = 5.8 s. Panel (b) indicates the thermal plasma profiles and panel (c) the EP density profiles used in the analysis of the AE stability with respect to the NBI deposition region. The iota profile assumed in the study is obtained from the MSE data shown in Fig. 4. Figure 13 shows the growth rate and frequency of AEs calculated in the simulation for t = 5.8 if the NBI is deposited on-axis with r_peak = 0.1 [panels (a) and (b)], off-axis with r_peak = 0.35 [panels (c) and (d)], or strongly off-axis with r_peak = 0.6 [panels (e) and (f)]. The simulations with on-axis NBI injection show the largest AE growth rate, 2 times larger compared to the simulation with r_peak = 0.35 and 4 times larger compared to r_peak = 0.6 cases. In addition, the EP β threshold is 0.005 if r_peak = 0.1, increasing to 0.01 if r_peak = 0.6. The dominant mode is a $1/1−1/2$ TAE with a frequency around 100 kHz that moves radially outward as the NBI is deposited more off-axis [see Fig. 14(a) and 14(b)]. Consequently, the configuration with off-axis NBI heating has an improved AE stability with respect to the on-axis case because the EP β threshold is higher. In addition, if the EP β threshold is exceeded, the growth rate of the unstable $1/1−1/2$ TAE is lower in the off-axis case. That means, a stronger NBI injection power can be applied to heat the plasma before the $1/1−1/2$ TAE is triggered. Regarding the EP energy, EP populations with T_f > 120 keV destabilize the AEs with the largest growth rate. Hence, the application of an NBI with a higher voltage and the generation of EP populations with birth energies above 180 keV may reduce the EP β threshold and the unstable AEs would show a larger growth rate.

This section is dedicated to further analyze the discharge 167 800 performed controlling the NBCD and thermal plasma density to improve the stability of PGDM and AEs during the shot. The effect of the reduced AE and PGDM activity on the thermal plasma confinement and plasma heating efficiency is analyzed with respect to the energy contained in the plasma, confinement time as well as neutron and EP fluxes.

Figure 12(a) shows a rather large MP signal at frequencies below 10 kHz during the ctr-NBCD phase, from t = 3.3 to 5.3 s, indicating the destabilization of PGDM. At the beginning of the co-NBCD phase, between t = 5.3 and 5.5 s, the PGDM stability improves even though the plasma β increases [see Fig. 3(a)]. On the other hand, AEs in the frequency range of 50 to 70 kHz are destabilized due to the thermal plasma density decrease as well as the enhanced NBI power injection by the NBIs 1 and 3. From t = 5.5 s, PGDM are stable due to the up-shift of the iota profile (1/2 rational surface is non-resonant) as well as the decrease in the thermal plasma β. The thermal plasma density decreases by $50%$ from t = 5.5 to 5.8 s leading to a further destabilization of the AEs. Nevertheless, the up-shift of the iota profile has a stabilizing effect on the AEs, reducing the AE activity from t = 6.0 s, almost stable once the co-NBCD intensity reaches 20 kA.

Figure 15 shows the continuum gaps at t = 5.0 (black), 5.5 (red), and 5.8 s (red) in the discharge 167 800. There is an increase in the gaps frequency range between t = 5.0 and 5.8 s due to the combined effect of the thermal ion density decrease, thermal electron temperature increment (particular the upper frequency limit of the BAE gap), and the iota profile up-shift by the co-NBCD. The continuum gaps at t = 5.0 s are narrow reason why no AE activity is observed in the experiment. On the other hand, at t = 5.5 s there is a TAE gap located in the middle plasma region between the frequency range $40−50$ kHz. At t = 5.8 s, a wide TAE gap covers all the minor radius in the frequency range above 80 kHz. The up-shift of the TAE gap frequency range along the discharge is consistent with the AE frequency up sweeping observed in the experiment.

In summary, the combination of co-NBCD and plasma density control leads to the stabilization of PGDM and a reduction of the AEs activity in the discharge 167 800, improving the discharge performance based on the increment of the neutron fluxes. Nevertheless, EP losses induced by AEs cause a reduction of the NBI heating efficiency and a decrease in the energy contained in the plasma; thus, the discharge improved performance is transitory.

The NBCD modifies the iota profile and the continuum gaps in LHD experiments; thus, the MHD stability thresholds change. The generation of a large ctr-NBCD causes a down-shift of the iota profile leading to configurations with slender gaps in frequency range and radial location, particularly at the inner plasma region. The AE stability improves due to an enhancement of the continuum and FLR damping effects that is reflected in an improved plasma heating performance by the NBIs. Nevertheless, a TAE gap appears in the middle-outer plasma region resulting in unstable TAEs if the NBI is deposited off-axis. Ctr-NBCD configurations also show a deterioration of the plasma confinement due to the destabilization of 1/2 PGDM in the plasma core, induced by the 1/2 rational surface that resonates in the inner plasma region where the stabilizing effect of the magnetic shear is weak. On the other hand, co-NBCD causes an up-shift of the iota profile leading to configurations with wider continuum gaps and a TAE gap located at higher-frequency ranges compared to ctr-NBCD configurations. The AE stability improves due to an enhancement of the EP FLR damping effect because high-frequency AEs have a more localized and narrow eigenfunction, thus the eigenfunction width is closer to the EP Larmor radius and the energy of the AE is more efficiently transferred. If the EP β threshold is exceeded, TAEs covering a large fraction of the normalized minor radius can be triggered. Nevertheless, previous studies have shown that the EP β threshold for these configurations is rather high and the AE growth rate decreases as the co-NBCD intensity enhances.²³ Co-NBCD configurations causes the stabilization of 1/2 PGDM if the minima of the iota profile is above 0.5 and 1/2 rational surface is non-resonant.

MSE diagnostic shows 1/2 PGDM is stabilized once the 1/2 rational goes non-resonant due to the up-shift of the iota profile above 0.5 if the co-NBCD is large enough, consistent with previous experimental observations.⁹⁶ It must be reminded that uncertainty of MSE measurements may require the verification of the analysis performed once the diagnostic accuracy improves. The down-shift of the iota profile leads to a continuum configuration with a slender gap in radial width and frequency range, closing the TAE gap between the inner-middle plasma regions if the iota profile minima decreases below 0.2. Consequently, the TAEs identified in the simulations show a growth rate $25%$ lower compared to configurations with iota profile minima of 0.3 and 0.1.

The set of experiments performed to analyze the effect of the thermal plasma density on the MHD stability identifies opposite optimization trends for PGDM and AE. Low-density discharges show stable PGDM because the averaged thermal β of the plasma is below the threshold for the PGDM destabilization, i.e., the plasma pressure gradient is not large enough to overcome the stabilizing effect of the magnetic shear and thermal ion FLR damping. On the other hand, AEs are unstable in discharges with low density because low-density configurations show wide continuum gaps, particularly a TAE gap between the inner-middle plasma regions that causes the destabilization of TAEs for rather low EP β values. The experiments show unstable AEs if the averaged thermal plasma density is lower than $2×1019$ m⁻³ although the PGDM are stable or marginally unstable if the averaged thermal plasma β is below 0.65. Configurations with an averaged plasma density between $1.2−2.1×1019$ m⁻³ indicate the co-existence of PGDM and AEs, i.e., plasma confinement and heating are handicapped by MHD instabilities. Consequently, such configurations must be avoided to improve the LHD performance. It should be noted that configurations with averaged plasma density between $1.7−2.7×1019$ m⁻³ show the destabilization of modes in the frequency range between $20−40$ kHz, identified in the simulations as PGDM coupled with shear Alfvénic waves triggered in the inner-middle plasma region.

High-density LHD operations above $3×1019$ m^{− 3} lead to the stabilization of AEs although a strong destabilization of PGDM. In particular, even though 1/2 PGDM is stable if a 1/2 rational surface goes non-resonant, peripheral modes as 1/1 PGDM are strongly destabilized in high thermal β discharges, above $2%$⁠. Consequently, another actuators should be applied to avoid or minimize the effect of PGDM on the high thermal β discharges confinement, for example, applying an off-axis injection of ECH or ICRH to reduce the plasma resistivity, the generation of ECCD to locally increase the stabilizing effect of the magnetic shear or the application of resonant magnetic perturbation (RMP) to avoid mode locking.^97–100 

The AE stability in LHD low-density plasma can be also optimized with respect to the NBI operational regime. The theoretical analysis performed for the discharge 167 800 at t = 5.8 s indicate off-axis NBI injection has an EP β threshold of 0.01, two times larger compared to configurations with on-axis NBI heating. Thus, the $1/1−1/2$ TAE identified in the simulations remains stable for a NBI injection power two times larger in operation scenarios with off-axis heating compared to the on-axis case. The analysis of the resonance intensity with respect to the EP energy indicates EP populations with energies above 120 keV cause the destabilization of the AEs with the largest growth rates. This is explained by an enhancement of the EP resonance as the EP energy increases and the ratio between EP and Alfvén velocity is closer to unity. Consequently, the application of an NBI with a higher voltage in the LHD plasma may generate EP populations with birth energies above 180 keV, decreasing the EP β threshold for AE destabilization leading to higher growth rates and reducing the NBI heating efficiency. It should be noted that the voltage of LHD NBIs cannot be modified so any optimization trend linked to the EP birth energy cannot be explored. Likewise, the deposition region of the NBI cannot be modified by changing the NBI tilt because the injectors are fixed. Nevertheless, the AE stability for an off-axis NBI injection can be analyzed by displacing the vacuum magnetic axis outward (LHD outwards shifted configurations), in scenarios with high thermal β and a large Shafranov shift as well as operations with large plasma density and poor beam penetration. It should be noted that the stiffness of the EP density profile was not included in the parametric study, but it has an important effect on the AE stability. Minimizing the EP density gradient increases the EP β threshold for AE destabilization. This is the case of NBI operational regimes leading to a large radial spreading of the beam deposition, avoiding a strongly localized beam deposition and an EP radial density profile with large gradients.

It must be recalled the numerical analysis performed is based on the linear stability of PGDM and AE. Nevertheless, the saturation phase of PGDM and AE must be also studied to fully understand the consequences of these instabilities on the device performance. Examples of adverse events observed during the saturation phase of PGDM in LHD inward-shifted discharges are the internal collapses induced by 1/2⁴⁰ or 1/1 PGDM,^101,102 sawtooth-like events caused by 1/2 or 1/3 PGDM,¹⁰³ and tongue-shaped deformations induced by 1/1 PGDM.^104,105 Outward shifted configurations also show core density collapses induced by unstable ballooning modes.^94,106,107 Such adverse events cause a partial de-confinement of the plasma although not necessarily the termination of the discharge. Likewise, MHD bursts are observed during the saturation phase of TAEs¹⁰⁸ and EIC burst during the saturation phase of the energetic-ion-driven resistive interchange mode,⁴⁷ leading to large EP losses before thermalization and inefficient plasma heating.^109–112 

Nonlinear simulations indicate the radial overlapping between magnetic islands induced by different PGDM may generate wide magnetic field stochastic regions, leading to the partial plasma de-confinement during internal collapse^81,83 and sawtooth-like events.^80,82 Nonlinear simulations also show MHD bursts are caused by the overlapping of n = 1 to 3 TAEs⁸⁴ and the EIC burst by the overlapping of the 1/1 EIC with 3/4 and 2/3 EPMs.⁶⁶ Consequently, fusion devices may have a large decrease in the plasma confinement and heating efficiency during the saturation phase of PGDM and AE if there is a transition from the soft to the hard MHD limit, i.e., a regime showing instabilities with an amplitude large enough to induce global plasma relaxations due to resonance overlapping. Thus, future nuclear fusion reactors may operate in the soft MHD limit (no overlapping between resonances).

The effect of the magnetic field intensity on the stability of PGDM and AE was not explored in the experiments, and it is not included in the analysis. Nevertheless, the magnetic field intensity is another important actuator for the identification of optimized scenarios. In particular, high-field devices show a larger Alfvén velocity leading to a lower ratio between EP and Alfvén velocities. Consequently, the design of an NBI that operates using a super-Alfvénic beam to minimize the destabilizing effect of high-energy EPs may require a large voltage.¹¹³ In addition, the target operation density must be also included in the NBI design because high-density operation leads to a decrease in the Alfvén velocity, although the scaling of the Alfvén velocity with the plasma density is weaker with respect to the magnetic field intensity. Also, device operation at a higher plasma density to weaken the EP resonance and reduce the AE activity is a strategy that can be easily applied in high-field Stellarator, although high-field Tokamaks are limited by the Greenwald density limit.¹¹⁴ 

A set of experiments during the 22th and 23th LHD campaigns were performed to analyze the stability of AEs and PGDM in inward-shifted LHD discharges. For the first time, the stability trends of AE and PGDM are analyzed at the same time to identify operation regimes showing an improved thermal plasma confinement and plasma heating performance. The actuators in the analysis are the NBI current drive, thermal plasma density, and NBI operational regime. MSE diagnostic provides information of the iota profile evolution during the co-NBCD phases of the discharges. AE and PGDM stability trends are studied using the FAR3d code with respect to the magnetic configuration, thermal plasma density, and NBI operational regime (NBI power injection, voltage and injection radial location).

Experimental data and simulation results indicate co-NBCD causes the stabilization of the 1/2 PGDM because the up-shift of the iota profile induces a transition of the 1/2 rational surface from resonant to non-resonant. On the other hand, ctr-NBCD leads to an improved AE stability due to a down-shift of the iota profile and a reduction of the frequency and radial width of the Alfvén gaps.

An increase in the thermal plasma density leads to the AEs stabilization because the EP resonance is weakened, the frequency range of the Alfvén gaps is slender, and the EP slowing-down time is smaller (lower EP β). The payback is a strong destabilization of the PGDMs.

The modification of the iota profile is more efficient to improve the AE stability compared to an increment of the thermal plasma density. The down-shift of the iota profile reduces frequency and radial width of the Alfvén gaps although, increasing the thermal plasma density, the gap radial width is weakly affected.

The NBI operational regime can be optimized to minimize the AE activity if the beam is deposited off-axis, leading to a higher EP β threshold to destabilize AEs.

The authors would like to thank the LHD technical staff for their contributions in the operation and maintenance of LHD. This work was supported by the Comunidad de Madrid under the Project No. 2019-T1/AMB-13648, U.S. DOE under Grant No. DE-FG02-04ER54742, and NIFS07KLPH004. The data supporting the findings of this study are available in the LHD experiment data repository at https://doi.org/10.57451/lhd.analyzed-data.

The authors have no conflicts to disclose.

Jacobo Varela: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Donald A. Spong: Investigation (supporting); Software (equal); Writing – review & editing (supporting). Luis Garcia: Investigation (supporting); Software (equal); Writing – review & editing (supporting). Yashika Ghai: Investigation (supporting); Software (supporting); Writing – review & editing (supporting). Juan Ortiz: Investigation (supporting). Kenichi Nagaoka: Data curation (supporting); Investigation (supporting); Writing – review & editing (supporting). Yuki Takemura: Data curation (supporting); Investigation (supporting); Writing – review & editing (supporting). Kiyomasa Watanabe: Data curation (supporting); Investigation (supporting); Writing – review & editing (supporting). Katsumi Ida: Data curation (supporting); Investigation (supporting); Writing – review & editing (supporting). Mikiro Yoshinuma: Data curation (supporting); Investigation (supporting). Kazunobu Nagasaki: Data curation (supporting); Investigation (supporting); Writing – review & editing (supporting). Alvaro Cappa: Data curation (supporting); Investigation (supporting); Writing – review & editing (supporting). Sergei E. Sharapov: Data curation (supporting); Investigation (supporting); Writing – review & editing (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.