Subdominant modes and optimization trends of DIII-D reverse magnetic shear configurations

TL;DR

This study uses advanced simulations to analyze MHD and Alfven Eigenmode activity in DIII-D reverse magnetic shear configurations, revealing how different operational regimes influence plasma stability and mode destabilization.

Contribution

It introduces a comprehensive 3D reduced MHD model coupled with energetic particle dynamics to identify conditions minimizing MHD and AE activity in reverse shear plasmas.

Findings

Reverse shear regions destabilize ideal MHD modes with higher growth rates.

Off-axis NBI injection improves AE stability compared to on-axis injection.

Adjusting q profile and shear near the axis can enhance MHD stability.

Abstract

Alfven Eigenmodes and magneto-hydrodynamic modes are destabilized in DIII-D reverse magnetic shear configurations and may limit the performance of the device. We use the reduced MHD equations in a full 3D system, coupled with equations of density and parallel velocity moments for the energetic particles (with gyro-fluid closures) as well as the geodesic acoustic wave dynamics. The aim of the study consists in finding ways to avoid or minimize MHD and AE activity for different magnetic field configurations and neutral beam injection operational regimes. The simulations show at the beginning of the discharge, before the reverse shear region is formed, a plasma that is AE unstable and marginally MHD stable. As soon as the reverse shear region appears, ideal MHD modes are destabilized with a larger growth rate than the AEs. Both MHD modes and AEs coexist during the discharge, although the MHD modes are more unstable as the reverse shear region deepens. The simulations indicate the destabilization of Beta induced AE, Toroidal AE, Elliptical AE and Reverse Shear AE at different phases of the discharges. A further analysis of the NBI operational regime indicates that the AE stability can be improved if the NBI injection is off axis, because on-axis injection leads to AEs with larger growth rate and frequency. In addition, decreasing the beam energy or increasing the NBI relative density leads to AEs with larger growth rate and frequency, so an NBI operation in the weakly resonant regime requires higher beam energies than in the experiment. The MHD linear stability can be also improved if the reverse shear region and the q profile near the magnetic axis are in between the rational surfaces q=2 and q=1, particularly if there is a region in the core with negative shear, avoiding a flat q profile near the magnetic axis.

Figures

Tables

Table 1. Table 1: Plasma parameters (at the magnetic axis) at different times along the discharges. First column is the thermal ion density, second column is the thermal ion temperature, third column is the thermal β 𝛽 \beta , forth column is the energetic particle density, fifth column is the energetic particle temperature, sixth column the energetic particle β 𝛽 \beta and seventh column the normalized cyclotron frequency.

Shot $164841$ (Low frequency AEs)

Table 2. Table 2: Dynamic and equilibrium toroidal (n) and poloidal (m) modes.

Dynamic (m) modes

Equations

$$\mathbf{v}=\sqrt{g}{R}_0\nabla \zeta \times \nabla \Phi ,\mathbf{B}=R_0\nabla \zeta \times \nabla \psi ,$$

$$\frac{\partial \tilde{\psi }}{\partial t}=\sqrt{g}B{\nabla }_{\Vert }\Phi +\eta {\epsilon }^{2}J{\tilde{J}}^{\zeta }$$

$$\frac{\partial \tilde{U}}{\partial t}=-ϵv_{\zeta ,eq}\frac{\partial U}{\partial \zeta }$$

$$+S^2\left[\sqrt{g}B{\nabla }_{\Vert }{J}^{\zeta }-\frac{{\beta }_0}{2{\epsilon }^2}\sqrt{g}{\left(\nabla \sqrt{g}\times \nabla \tilde{p}\right)}^{\zeta }\right]$$

$$-S^2\left[\frac{{\beta }_f}{2{\epsilon }^2}\sqrt{g}{\left(\nabla \sqrt{g}\times \nabla {\tilde{n}}_f\right)}^{\zeta }\right]$$

$$\frac{\partial \tilde{p}}{\partial t}=-ϵv_{\zeta ,eq}\frac{\partial p}{\partial \zeta }+\frac{d{p}_{eq}}{d\rho }\frac{1}{\rho }\frac{\partial \tilde{\Phi }}{\partial \theta }$$

$$+\Gamma p_{eq}\left[\sqrt{g}{\left(\nabla \sqrt{g}\times \nabla \tilde{\Phi }\right)}^{\zeta }-{\nabla }_{\Vert }{v}_{\Vert th}\right]$$

$$\frac{\partial {\tilde{v}}_{\Vert th}}{\partial t}=-ϵv_{\zeta ,eq}\frac{\partial v_{||th}}{\partial \zeta }-\frac{S^2{\beta }_0}{n_{0,th}}{\nabla }_{\Vert }p$$

$$\frac{\partial {\tilde{n}}_f}{\partial t}=-ϵv_{\zeta ,eq}\frac{\partial n_f}{\partial \zeta }-\frac{S{v}_{th,f}^2}{{\omega }_{cy}}{\Omega }_d({\tilde{n}}_f)-S{n}_{f0}{\nabla }_{\Vert }{v}_{\Vert f}$$

$$-{\epsilon }^2{n}_{f0}{\Omega }_d(\tilde{\Phi })+{\epsilon }^2{n}_{f0}{\Omega }_{\ast }(\tilde{\Phi })$$

$$\frac{\partial {\tilde{v}}_{\Vert f}}{\partial t}=-ϵv_{\zeta ,eq}\frac{\partial v_{||f}}{\partial \zeta }-\frac{S{v}_{th,f}^2}{{\omega }_{cy}}{\Omega }_d({\tilde{v}}_{\Vert f})$$

$$-{\left(\frac{\pi }{2}\right)}^{1/2}S{v}_{th,f}\left|{\nabla }_{\Vert }{v}_{\Vert f}\right|$$

$$-\frac{S{v}_{th,f}^2}{n_{f0}}{\nabla }_{\Vert }{n}_f+S{\epsilon }^2{v}_{th,f}^2{\Omega }_{\ast }(\tilde{\psi })$$

$$J^{\zeta }=\frac{1}{\rho }\frac{\partial }{\partial \rho }\left(-\frac{g_{\rho \theta }}{\sqrt{g}}\frac{\partial \psi }{\partial \theta }+\rho \frac{g_{\theta \theta }}{\sqrt{g}}\frac{\partial \psi }{\partial \rho }\right)$$

$$-\frac{1}{\rho }\frac{\partial }{\partial \theta }\left(\frac{g_{\rho \rho }}{\sqrt{g}}\frac{1}{\rho }\frac{\partial \psi }{\partial \theta }+\rho \frac{g_{\rho \theta }}{\sqrt{g}}\frac{\partial \psi }{\partial \rho }\right)$$

$${\Omega }_d=\frac{1}{2B^4\sqrt{g}}\left[\left(\frac{I}{\rho }\frac{\partial B^2}{\partial \zeta }-J\frac{1}{\rho }\frac{\partial B^2}{\partial \theta }\right)\frac{\partial }{\partial \rho }\right]$$

$$-\frac{1}{2B^4\sqrt{g}}\left[\left(\rho {\beta }_{\ast }\frac{\partial B^2}{\partial \zeta }-J\frac{\partial B^2}{\partial \rho }\right)\frac{1}{\rho }\frac{\partial }{\partial \theta }\right]$$

$$+\frac{1}{2B^4\sqrt{g}}\left[\left(\rho {\beta }_{\ast }\frac{1}{\rho }\frac{\partial B^2}{\partial \theta }-\frac{I}{\rho }\frac{\partial B^2}{\partial \rho }\right)\frac{\partial }{\partial \zeta }\right]$$

$${\Omega }_{\ast }=\frac{1}{B^2\sqrt{g}}\frac{1}{n_{f0}}\frac{d{n}_{f0}}{d\rho }\left(\frac{I}{\rho }\frac{\partial }{\partial \zeta }-J\frac{1}{\rho }\frac{\partial }{\partial \theta }\right)$$

$${\nabla }_{\Vert }f=\frac{1}{B\sqrt{g}}\left(\frac{\partial \tilde{f}}{\partial \zeta }+to0.0pt-\text{hss}\iota \frac{\partial \tilde{f}}{\partial \theta }-\frac{\partial f_{eq}}{\partial \rho }\frac{1}{\rho }\frac{\partial \tilde{\psi }}{\partial \theta }+\frac{1}{\rho }\frac{\partial f_{eq}}{\partial \theta }\frac{\partial \tilde{\psi }}{\partial \rho }\right)$$

$$\sqrt{g}{\left(\nabla \sqrt{g}\times \nabla \tilde{f}\right)}^{\zeta }=\frac{\partial \sqrt{g}}{\partial \rho }\frac{1}{\rho }\frac{\partial \tilde{f}}{\partial \theta }-\frac{1}{\rho }\frac{\partial \sqrt{g}}{\partial \theta }\frac{\partial \tilde{f}}{\partial \rho }$$

$$\frac{1}{\sqrt{g}}=\frac{B^2}{{\epsilon }^2(J+to0.0pt-\text{hss}\iota I)}$$

$$n_f(r)=\frac{(0.5(1+\tanh (r_{flat}\ast (r_{peak}-r))+0.02)}{(0.5(1+\tanh (r_{flat}\ast r_{peak}))+0.02)}$$

Full text

Keywords: Tokamak, DIII-D, Pedestal, MHD, AE, energetic particles

This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05- 00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

\ioptwocol

1 Introduction

Reverse magnetic shear discharges are a possible operational scenario of nuclear fusion reactors because these configurations have good MHD stability properties and energy confinement at high $\beta$ [1, 2, 3, 4, 6, 7], although Alfvén Eigenmode stability and energetic particle confinement are still an open issue [8, 9, 10, 11, 12, 13, 14]. MHD stability of interchange and ballooning modes is improved in high $\beta$ discharges by combining reverse shear operation (the plasma is in the second stable regime [15, 16]) and fast toroidal rotation [17]. In addition, the energy confinement is improved in reverse shear configurations due to the stabilization of drift-type micro-instabilities [18]. Also, significant bootstrap current can be generated, lowering the current drive requirement for steady-state operation [19, 20].

Reverse magnetic shear discharges are extensively studied in the DIII-D device [21, 22, 23] showing a large fraction of bootstrap current [24, 25, 26] and improved MHD stability [27, 28, 29]. These configurations have been selected as a base line scenario for ITER and DEMO [30, 31, 32], although an enhancement of the energetic particle transport by unstable Alfven Eigenmodes was measured [33, 34, 35].

The transport of energetic neutral beam ions, particles heated using ion cyclotron resonance heating (ICRF) and fusion produced alpha particles can be enhanced by energetic particles driven instabilities [36, 37, 38]. Experiments in TFTR, JET and DIII-D tokamaks as well as LHD and W7-AS stellarators identified a decrease in the device performance caused by the destabilization of Alfven Eigenmodes [39, 40, 41, 42, 43, 44], leading to an increase of particle and diffusive loses if the transit, bounce or drift frequency of the energetic particles resonate with AE mode frequencies.

Super-Alfvénic alpha particles and energetic particles destabilize Alfvén Eigenmodes (AE) in the spectral gaps of the shear Alfvén continua [45, 46]; these have been seen in several DIII-D configurations [47, 48, 49, 50]. Periodic variations of the Alfvén speed destabilize different AE families, for example toroidicity induced Alfvén Eigenmodes (TAE) coupling $m$ with $m+1$ modes (m is the poloidal mode) [51, 52], ellipticity induced Alfvén Eigenmodes (EAE) coupling $m$ with $m+2$ modes [53, 54] or noncircularity induced Alfvén Eigenmodes (NAE) coupling $m$ with $m+3$ or higher [55, 56]. Other example of AE are the beta induced Alfvén Eigenmodes (BAE) driven by compressibility effects [57], Reversed-shear Alfvén Eigenmodes (RSAE) due to local maxima/minima in the safety factor $q$ profile [33, 58] and Global Alfvén Eigenmodes (GAE) observed in the minimum of the Alfvén continua [59, 60]. Furthermore, current driven instabilities as internal kinks [61, 62] or pressure gradient driven instabilities as interchange and ballooning modes [63] can be kinetically destabilized. textcolorredPlease, see the Appendix for further details about the instability characteristics and eigenfunction structure, as well as a comparison with the electron temperature fluctuations (ECE) data measured during the discharges.

DIII-D has eight neutral beam injectors (NBI), six sources injected in the midplane (on axis) and 2 injected downwards at an angle (off axis) to heat, fuel and drive current in the plasma. Six sources are co-injected (direction of the plasma current) with two tilted sources and 2 source are counter-injected (opposite to the plasma current). The NBI injects deuterium with a beam energy of $80$ keV ($2.25$ MW source) in a deuterium plasma. Strong NBI heating leads to the destabilization of AE as TAE [64], EAE [65], NAE [66], BAE [67], RSAE [68] and GAE [69]. DIII-D experiments measured a larger transport and enhanced energetic particle losses linked to unstable AE [70, 71, 72].

The aim of the present study is to identify optimization pathways to improve MHD and AE linear stability of DIII-D reverse magnetic shear discharges, comparing simulation results and experimental observations. We analyze the effect of the magnetic field configuration and NBI operation regime (beam energy, injection intensity and deposition profile) on MHD and AE stability. We identify dominant and sub-dominant MHD and AE modes at different phases in a selection of three reverse magnetic shear DIII-D discharges. In addition, we study the overlapping between low frequency AE (frequency smaller than $40$ kHz) and MHD modes, as well as the effect of the energetic particle forcing on the MHD modes growth rate and frequency.

This analysis is performed using the FAR3D code [73, 74, 75], with extensions to include the moment equations of the energetic ion density and parallel velocity [76, 77] allowing treatment of linear wave-ion resonances. The numerical code solves the reduced non-linear resistive MHD equations adding the Landau damping/growth (wave-particle resonance effects) and geodesic acoustic waves (parallel momentum response of the thermal plasma) [58]. The Landau closure relations (Landau-closure models can be calibrated by more complete kinetic models [58]). The simulations are based on an equilibria calculated by the VMEC code [78].

This paper is organized as follows. The model equations, numerical scheme and equilibrium properties are described in section 2. The study of dominant and subdominant modes at different phases of reverse magnetic shear discharges is shown in section 3. The optimization trends to minimize AE and MHD activity modifying the NBI operational regime and the magnetic field configuration are analyzed in section 4. Studies of the interaction between MHD modes and low frequency AE is presented in section 5. Finally, the conclusions of this paper are presented in section 6.

2 Equations and numerical scheme

From the full MHD equations we can derive a reduced set of three equations that follow the evolution of the pressure $p$, the stream function proportional to the electrostatic potential $\Phi$ and the perturbation of the poloidal flux $\psi$ while retaining the toroidal angle variation $\zeta$ (exact three-dimensional equilibrium) [79]. The equations for the background or thermal plasma evolution are complemented by moments of the energetic ions to add the effect of the wave-particle interaction [80], namely the energetic particle density ($n_{f}$) and velocity moments parallel to the magnetic field lines ($v_{||f}$). The coefficients of the closure relation are selected to match a two-pole approximation of the plasma dispersion function.

The reduced equations are based on the following assumptions: high aspect ratio, medium $\beta$ (of the order of the inverse aspect ratio $\varepsilon=a/R_{0}$), small variation of the fields and small resistivity. The magnetic field and plasma velocity perturbations are defined as:

[TABLE]

The equations, in dimensionless form, are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here, $U=\sqrt{g}\left[{\nabla\times\left({\rho_{m}\sqrt{g}{\bf{v}}}\right)}\right]^{\zeta}$ is the vorticity and $\rho_{m}$ the ion and electron mass density. The toroidal current density $J^{\zeta}$ is defined as:

[TABLE]

$v_{\zeta,eq}$ is the equilibrium toroidal rotation and $v_{||th}$ is the parallel velocity of the thermal particles. $n_{f}$ is normalized to the density at the magnetic axis $n_{f_{0}}$, $\Phi$ to $a^{2}B_{0}/\tau_{R}$ and $\Psi$ to $a^{2}B_{0}$. All lengths are normalized to a generalized minor radius $a$; the resistivity to $\eta_{0}$ (its value at the magnetic axis); the time to the resistive time $\tau_{R}=a^{2}\mu_{0}/\eta_{0}$; the magnetic field to $B_{0}$ (the averaged value at the magnetic axis); and the pressure to its equilibrium value at the magnetic axis. The Lundquist number $S$ is the ratio of the resistive time to the Alfvén time $\tau_{A0}=R_{0}(\mu_{0}\rho_{m})^{1/2}/B_{0}$. $\hbox to0.0pt{-\hss}\iota$ is the rotational transform, $v_{th,f}=\sqrt{T_{f}/m_{f}}$ the energetic particle thermal velocity normalized to the Alfvén velocity in the magnetic axis $v_{A0}$ and $\omega_{cy}$ the energetic particle cyclotron frequency times $\tau_{A0}$. $q_{f}$ is the charge, $T_{f}$ the temperature and $m_{f}$ the mass of the energetic particles. The $\Omega$ operators are defined as:

[TABLE]

[TABLE]

Here the $\Omega_{d}$ operator is constructed to model the average drift velocity of a passing particle and $\Omega_{*}$ models its diamagnetic drift frequency. We also define the parallel gradient and curvature operators:

[TABLE]

[TABLE]

with the Jacobian of the transformation:

[TABLE]

The geodesic compressibility in the frequency range of the geodesic acoustic mode (GAM) is included in the model by the parallel momentum response of the thermal plasma in the equations 2 and 5 [81, 82].

The equations are written using the equilibrium flux coordinates $(\rho,\theta,\zeta)$ with $\rho$ the generalized radial coordinate proportional to the square root of the toroidal flux function (normalized to one at the edge) and $\theta$ the poloidal angle. We use the Boozer formulation for the flux coordinates [83] with $\sqrt{g}$ the Jacobian of the coordinate transformation. The functions have two components, equilibrium and perturbation, represented as: $A=A_{eq}+\tilde{A}$.

The FAR3D code uses finite differences in the radial direction and Fourier expansions in the two angular variables. The numerical scheme is semi-implicit in the linear terms. The nonlinear version uses a two semi-step method to ensure $(\Delta t)^{2}$ accuracy.

Many of the results in this paper have been obtained utilizing an eigenvalue solver, which is a recently added feature to the FAR3D model. In the linear regime, the time variation of all quantities will asymptote to $e^{-i\omega t}$; in this limit the equations of the gyrofluid model can be cast into the form of a generalized eigenvalue problem: $\omega Ax=Bx$. This problem can be solved using an extension of the Jacobi-Davidson algorithm that was applied in [84] for stable shear Alfvén eigenmodes to complex eigenvalues. For the complex case the user provides a target value for the frequency and growth rate and the algorithm returns a number of eigenvalues that are near the target in the complex plane. This method is especially useful in searching for higher frequency Alfvénic instabilities when lower frequency MHD instabilities are present with higher growth rates. If the initial value approach was used, only the faster growing MHD modes would be visible. The identification of subdominant and damped modes is also becoming an important tool [85] in understanding the nonlinear dynamics of plasma instabilities.

The model was validated on previous studies of AE activity in LHD [86, 87], TJ-II [88, 89, 90] and DIII-D [91, 92] indicating a reasonable agreement with the experimental observations, identifying the main instability properties of the discharges, that is to say, the instabilities associated with the dominant modes and most important sub-dominant modes. It should be noted that the destabilizing effect of the thermal ion temperature gradients are not included in the model, so the Alfvenic ion temperature gradient (AITG) instability is out of the scope of the present analysis [93, 94]. The only thermal plasma drive included here is through the total thermal plasma pressure gradient. However, this only provides a drive for MHD ballooning modes, not for ITG instabilities. A global gyrofluid model has been constructed for ITG modes [95] and, in principle, such drives could be included in this model, but this will remain a topic for future research.

The EP and thermal ion finite Larmor radius damping effects as well as the electron-ion Landau damping effects are not included in the study for simplicity, leading to simulations with larger growth rates as compared with the experimental observations. Nevertheless, the simulations reproduce the main instabilities observed in the experimental data. Also, an analysis based on calculations using the instability growth rates is more robust from the point of view of the optimization as discussed in section 4.

2.1 Equilibrium properties

The table 1 shows the shot numbers and time frames of the fixed boundary VMEC equilibrium used in the simulations, calculated from the DIII-D reconstruction of reverse magnetic shear discharges [78].

The shot $164841$ is a case with strong low frequency AE activity, shot $164842$ a discharge with central positive shear and shot $164922$ a discharge with central negative shear. The equilibrium reconstruction is done using the following experimental constrains: magnetic data, MSE data, kinetic pressure and edge density profile from NEO model. We use up-down symmetric equilibria, although deviations from the original single-null divertor equilibria are small. We analyze four different phases during the discharges: phase 0 before the formation of the reverse shear region, phase A during the destabilization of AE activity with steady frequencies, phase B during the destabilization of AE activity with non steady frequencies and phase C before the ramp down. Table 1 summarizes the main plasma parameters. The magnetic field at the magnetic axis is $1.87$ T, the averaged inverse aspect ratio is $\varepsilon=0.48$ and the average elongation is $\Delta=1.7$.

Figure 1 shows the Alfvén gaps of $n=4$ toroidal mode of shot $164842$ at the different discharge phases. There are four main Alfvén gaps: below $50$ kHz (BAE, BAAE and GAE gap), between $[50,175]$ kHz (TAEs gap), between $[175,350]$ kHz (EAE gap) and above $350$ kHz (NAE gap). The other shots analyzed in present study show similar gap distribution.

2.2 Simulation parameters

Table 2 shows the dynamic and equilibrium toroidal (n) and poloidal (m) modes for each model. The plasma edge ($\rho<0.75$) is not included in the analysis of M2 and M3 cases to avoid he excitation of artificial ideal MHD instabilities at the outer boundary. The mode number notation is $m/n$ consistent with the $q=m/n$ definition. The simulations have a uniform radial grid of 1000 points (200 points for the subdominant modes studies).

The MHD parities are broken by the kinetic closure moment equations (6) and (7) so both parities are included in the study. For all the dynamic variables $sin(m\theta+n\zeta)$ and $cos(m\theta+n\zeta)$ parities are considered in the calculation of the growth rate and real frequency. The convention is as follows: $n>0$ corresponds to the Fourier component $\cos(m\theta+n\zeta)$ and $n<0$ to $\sin(-m\theta-n\zeta)$ (pressure eigenfunction case). The magnetic Lundquist number is $S=5\cdot 10^{6}$ similar to the experimental value in the middle of the plasma.

In the simulations the NBI injection intensity is controlled by fast ion $\beta$ value ($\beta_{f}$), linked to the density ratio between energetic particles and bulk plasma $n_{f}(0)/n_{e}(0)$ at the magnetic axis, calculated by the code TRANSP without the effect of the anomalous beam ion transport. The beam energy and the resonance regime between thermal plasma and energetic particles is controlled by the ratio of the energetic particle thermal velocity to the Alfvén velocity at the magnetic axis ($v_{th,f}/v_{A0}$). The energetic particle distribution function is a Maxwellian.

The shots analyzed use co-injected NBI. These discharges are categorized as negative shear cases although during the shot $164922$ there are phases where the safety factor shows a positive shear near the magnetic axis, reason why we identify this discharge as a positive shear case [96].

3 Dominant and subdominant MHD modes and AE

In this section we analyze the growth rate ($\gamma$) and frequency ($f$) of the dominant and sub-dominant MHD and AE modes. We compare the simulation results with the instabilities observed in the experiments by magnetic and CO2 interferometry diagnostics. The instabilities with the largest growth rate that can constrain the device performance are identified in each phase of the discharge.

Low frequency AEs

Figure 2 shows the instabilities measured during shot $164841$. In this discharge both CO2 interferometer and magnetic data are available. The CO2 interferometry shows AE activity throughout the discharge: burst activity at M10, constant frequency AE at M1A, up-sweeping frequency AE at M1B and weaker steady frequency AE at M1C. In addition, the magnetic diagnostic shows $n=1$ to $4$ low frequency instabilities in the range of $5$ to $40$ kHz at different phases of the discharge. In figure 2 the colored stars show frequencies of the dominant modes in the simulations overlaid on the experimental diagnostic signals. As can be seen the simulations bracket many of the experimentally measured frequency lines. However, they miss some of the higher frequency lines for the M10, M1A, and M1B time slices. Also, in some cases modes are predicted at lower frequencies (M10 and M1B) than seen experimentally.

Figure 3 shows the plasma profiles and magnetic field configuration at different phases of the discharge. It should be noted that the energetic particle density gradient and the reverse shear region are both located between the inner and middle plasma, so RSAE are easily destabilized.

Figure 4 shows the growth rate and frequency of dominant and subdominant modes. At M10 phase, only n=1 modes are MHD unstable although all n modes are AE unstable: $n=1$ BAE with $f\approx 30$ kHz and $n=2$ to $n=6$ TAE with $f\approx 55$, $75$, $120$, $150$ and $160$ kHz. The AEs growth rate is between 2 to 3 times larger than that of the MHD modes. This result is consistent with the observation that the discharge shows strong AE activity in its early phase. At M1A phase, all modes are MHD unstable with growth rates 2 and 3 times larger than the unstable AE: $n=1$ and $n=2$ BAE with $f\approx 25$ kHz, $n=4$ BAE with $f\approx 45$ kHz, $n=3$ to $n=6$ TAE with $f\approx 75$, $105$, $125$ and $150$ kHz and $n=5$ and $6$ EAE with $f\approx 175$ and $195$ kHz. The strong $n=2$ instability (with the $n=4$ overtone) observed by the magnetic diagnostic is reproduced in the simulations as $n=2$ and $n=4$ BAE, as well as the $n=5$ and $n=6$ TAE/EAE, consistent with the range of instability frequencies measured by CO2 interferometer. The enhancement of the MHD activity is linked to the formation of the reverse shear region in the inner-middle plasma, as well as an increase of the AE growth rate by around $50\%$. At M1B phase, again the most unstable modes are MHD but the growth rate is slightly smaller compared to M1A phase, although the growth rate of the AEs is similar for the $n=1$ BAE with $f\approx 30$ kHz, the $n=2$ to $n=6$ TAE with $f\approx 60$, $95$, $125$, $145$ and $160$ kHz as well as $n=5$ and $6$ EAE with $f\approx 175$ and $190$ kHz, showing an increase in the AE frequency compared to M1A phase. The main characteristic of this discharge phase is the destabilization of RSAEs, resulting in frequencies in the range of the TAE gap and growth rates slightly larger than the TAEs, consistent with the up-sweeping frequency AEs observed in the experiment. The destabilization of the RSAE is caused by a slight increase of the EP density gradient and a decrease of the reverse shear region depth. At M1C phase, the most unstable modes are still MHD as well as the subdominant $n=2$ and $3$ BAE with $f\approx 25$ and $30$ kHz, similar to the magnetic diagnostic observations. In addition $n=1$ BAE with $f\approx 45$ kHz, $n=2$ to $n=6$ TAE with $f\approx 60$, $90$, $110$, $140$ and $155$ kHz as well as $n=5$ and $6$ EAE with $f\approx 175$ and $200$ kHz are unstable. It should be noted that the simulations identify instabilities in the frequency range of $f=[70,125]$ kHz not observed in the experimental data. The instabilities calculated in the simulations in that frequency range are TAE, RSAE and BAE destabilized between the middle and the outer plasma region where the Alfven gap is extended to higher/lower frequencies, as well as sub-dominant energetic particles modes (EPM) destabilized in the continuum. The reason these modes are unstable in the simulations but not observed in the experiment is because the model doesn’t include the EP finite Larmor radius or electron-ion Landau damping effects, that should lead to the stabilization of these modes. Future studies will include these damping effects to improve the simulation results. In addition, the safety factor and energetic particle profiles are not directly measured during the experiment, affecting the AE stability threshold obtained in the simulations. Newer fast ion Dα (FIDA) measurements in DIII-D show that the action of AE modes is to flatten the EP profiles, thus the EP density profiles could be less peaked than the configuration analyzed.

In summary, DIII-D performance in the early phase of the M1 discharge is limited by the AE activity driven in the inner plasma region before the reverse shear is formed. After the formation of the reverse shear region, MHD-likes modes are dominant and coexist with subdominant AEs. During the late phase of the discharge, RSAE are destabilized when the EP density gradient reaches it maximum value inside the reverse shear region. The optimization of the discharge stability based on these calculations would require modified phasing of the NBI power, particularly at the beginning of the discharge, and the amelioration of the MHD stability as soon as the reverse shear region is formed. The reason why this discharge has a bad MHD stability is because the inner plasma core resonates with the $q=3$ rational surface and the bottom of the reverse shear region with $q=2$ rational surface.

A more detailed discussion of how the AE families are identified by analyzing the individual eigenfunctions is given in the Appendix, as well as a comparison with the experimental data.

Central negative shear

Figure 5 shows the instabilities measured during shot $164842$. CO2 interferometry shows AE instabilities during all the discharge: burst activity at M20, constant frequency AE at M2A, up-sweeping frequency AE at M2B and down-sweeping frequency AE at M2C. The frequencies of the dominant modes in the simulations are consistent with the experimentally measured frequency lines, although some higher frequency lines are missing.

Figure 6 shows the plasma profiles at different phases of the discharge. The configuration is similar to the $164841$ shot although the $q$ profile near the magnetic axis is non constant and there is a negative shear region in the core.

Figure 7 shows the growth rate and frequency of dominant and subdominant modes at different phases of the discharge. The phases M20 and M2A show a similar AE stability than the M1 discharge, although the MHD stability is improved with only $n=1$ to $n=2$ unstable MHD modes, even if the reverse shear region was already formed at the M20 phase. On the other hand, $n=5$ and $n=6$ RSAE are destabilized during both phases. At phase M2B, $n=1$ to $n=3$ modes are MHD unstable as well as TAE/EAE/NAE of all toroidal modes in the range of frequencies $f=[60,700]$ kHz. The AE instabilities are stronger compared to the M1 discharge and the previous stages of the shot because there is a strong gradient of the EP density that covers the inner and middle plasma region, leading to the destabilization of high frequency AEs near the magnetic axis. At M2C phase, $n=1$ to $3$ MHD modes are unstable but the AE growth rates and frequencies are smaller compared to M2B phase because the EP density gradient only covers the inner plasma region.

Consequently, the MHD activity of shot $164842$ is smaller compared to $164841$ shot because the $q=3$ rational surface is non resonant in the inner plasma neither is $q=2$ near the bottom of the reverse shear region. On the other hand, AEs are destabilized from the first stage of the discharge, particularly the RSAE, because the energetic particle density gradient is located inside the reverse shear region. The discharge MHD and AE stability can be further improved increasing the magnetic shear in the inner plasma region, avoiding the destabilization of MHD modes and AE near the magnetic axis. In addition, unstable RSAE can be avoided if the energetic particle density gradient is not located near the bottom of the reverse shear region. Phase M2B shows the largest AEs growth rates and frequencies because the NBI beam energy and injection intensity are large enough to destabilize EAE and NAE modes of several toroidal families in the inner plasma region.

Central positive shear

Figure 8 shows the instabilities measured during shot $164922$. CO2 interferometry shows AE instabilities throughout the discharge: down-sweeping frequency AE at M30, constant frequency AE at M3A, up-sweeping frequency AE at M3B and high constant frequency AE ($f>250$ kHz) at M3C. Again, the frequencies of the dominant modes in the simulations are similar to the experimentally measured frequency lines, missing some higher frequency lines particularly in the M3B and M3C discharge phases.

Figure 9 shows the plasma profiles. The averaged thermal ion temperature and density are larger compared to shot $164842$. The energetic particle density gradient is narrower and located in the inner plasma region (between $0.2<r<0.4$), the $V_{th,f}/V_{A0}$ ratio is larger and there is a region of negative shear in the plasma core.

Figure 10 shows the growth rate and frequency of dominant and subdominant modes. At phase M30, the plasma is MHD stable but several AEs are destabilized, particularly $n=1-2$ TAEs and $n=6$ NAE with large growth rates and $f\approx 95$, $145$ and $375$ kHz destabilized near the magnetic axis where the q profile is almost flat. Again, this indicates the tendency for a weak magnetic shear to deteriorate the AE stability of the plasma. At phase M3A, $n=1$ and $2$ MHD modes are destabilized with growth rates 2 times larger compared to the AEs, especially the TAEs show the largest growth rates while the NAE are marginally unstable because the EP density gradient is weaker close to the magnetic axis. At phase M3B, $n=1$ to $n=4$ modes are MHD unstable, as well as several subdominant AE, in particular RSAE are destabilized in the middle plasma region. Phases M3B and M3C show similar MHD and AE stability properties. The range of frequencies observed by the magnetic measurements of low frequency AE and CO2 interferometer data are consistent with the simulation results.

These simulations indicate that the performance at the first stage of the discharge, while the reverse shear region is not yet formed, is limited by the TAE and EAE destabilized near the magnetic axis . During the rest of the discharge, $n=1$ to $n=4$ interchange modes are unstable in the bottom of the reverse shear region driven by $q=2$ rational surface. The magnetic shear near the inner plasma region is almost flat during all the discharge phases except at M3B phase, where TAE/RSAE are destabilized in the plasma core and in the bottom of the reverse shear region. In addition, subdominant TAE, EAE and NAE are destabilized in a range of frequencies up to $400$ kHz. BAE and low frequency TAE reproduced in the simulation are also consistent with the magnetic diagnostics, showing instabilities in the range of $f<60$ kHz during all the discharge phases.

In summary, after identifying the main dominant and sub-dominant modes that can impact the DIII-D heating efficiency at the different phases of the discharges, as well as some possible optimization trends to improve the plasma MHD and AE stability, in the next section we analyze the effect of the NBI operational regimen and magnetic field configuration on the plasma stability.

4 Optimization trends to improve the MHD and AE stability

In this section, based on the previous studies, we identify some optimization trends to improve the AE and MHD stability of the plasma by modifying the NBI operational regime and the magnetic field configuration.

4.1 NBI operational regime

We perform a parametric analysis varying the injection intensity ($\beta_{f}$), beam energy ($V_{th,f}/V_{A0}$) and the location of the EP density gradient ($r_{peak}$) to find the optimal operational regime of the NBI that minimizes AE activity. To study the effect of the EP density gradient location on the instabilities growth rate and frequency we define the next analytic expression for the EP density:

[TABLE]

where $r_{flat}$ parameter controls the slope flatness of the profile (fixed to $r_{flat}=7$). The study is performed at phase B of the M2 and M3 discharges.

Figure 11 shows the dependencies of growth rate and frequency for dominant and subdominant modes at different NBI injection intensities, beam energies and deposition regions for the M2B phase. Increasing $\beta_{f}$ from $0.01$ to $0.03$ (panels a to f) only leads to a transition between AE families (TAE $\rightarrow$ NAE, highlighted by a dotted orange arrow) for the $n=4$ to $6$ modes, although the critical $\beta_{f}$ to destabilize $n=1$ EAE and $n=2-3$ NAE is smaller than $0.01$, so there is only a partial optimization of the discharge AE stability if the NBI injection intensity is smaller than $\beta_{f}=0.01$. Increasing the beam energy ($V_{th,f}/V_{A0}$ ratio, panels g to m) leads to a transition between AE families for $n=1$ to $n=4$ modes: from $n=1$ TAE to $n=1$ BAE, $n=2$ EAE to $n=2$ TAE, $n=3$ EAE to $n=3$ TAE and $n=4$ NAE to $n=4$ TAE. On the other hand, no transition is observed for $n=5$ and $6$ EAE, only a decrease of the AE growth rate and frequency that can be explained as a weaker resonance between the EP and the thermal plasma. Changing the location of the EP density gradient from $r_{peak}=0.1$ (on-axis NBI injection) to $r_{peak}=0.3$ and $0.5$ (off-axis NBI injection) leads to the stabilization of $n=1$ TAE, $n=2$ EAE and $n=3$ to $6$ NAE (panels n to s), although $n=1$ BAE, $n=2-5$ TAE and $n=6$ to $6$ EAE are destabilized with lower growth rates, so the AE stability of the plasma is optimized if the NBI is deposited off-axis.

Figure 12 shows the same analysis for the M3B phase. Modifying $\beta_{f}$ (panel a to f) or the beam energy (panels g to m) do not lead to the stabilization of AEs in this range of parameters; the growth rate and frequency are similar, so no optimization trends are observed. On the other hand, off axis NBI heating leads to stabilization of EAE/NAE (panels n to s).

In summary, an improved NBI operational scenario requires higher beam energy and off-axis NBI heating. Decrease the NBI injection energy doesn’t lead to a significant improvement because the NBI operates above the critical $\beta_{f}$ value. It should be noted that there is experimental evidence of weaker AE activity if the beam energy decreases for the same injection intensity, as well as an enhancement of low frequency AEs. Consequently, an optimized NBI operational regime should not necessarily be linked to an arbitrary increase of the beam energy. In optimized NBI operational regimes, the EAE/NAE are stable and the growth rate/frequency of the BAE/TAE is smaller compared to the discharges analyzed. Previous studies also concluded that if the beam energy is high enough, the NBI operation entering into the non resonant regime where BAE/TAE modes are also stabilized [85, 88].

4.2 Magnetic field configuration

In this section, we analyze the effect of the magnetic field configuration with respect to MHD and AE stability at M2B and M3B phases. To accomplish this, the $q$ profile is displaced by $\Delta q=[-0.75,0.75]$ (Figure 13 case M3B) and we analyze the growth rate (panel a) and frequency (panels b and c) of the dominant $n=1$ to $6$ modes. $\Delta q>0$ displacements lead to larger growth rate and frequency for all MHD and AE instabilities compared to $\Delta q<0$ displacements. In addition, $\Delta q<0$ displacements lead to operation windows with improved AE and MHD stability (for example $\Delta q=-0.6$ compared to $\Delta q=-0.465$).

To confirm the AE and MHD stability improvement we analyze the dominant and subdominant modes growth rate and frequency if the q profile is displaced $\Delta q=-1$ in M2B case and $\Delta q=-0.6$ in M3B case, see Figure 13 panels d) and e), corresponding to operation windows with optimized stability properties. The improved M2B case shows weaker MHD activity with only an $n=1$ mode being unstable at lower growth rate. In addition, BAEs are marginally unstable and the dominant modes are $n=2-5$ TAEs and $n=6$ EAE. The upgraded version of M3B also shows weaker MHD activity for all modes (except $n=1$, slightly larger), although the growth rate and frequency of the AEs is similar to the original cases.

Consequently, the improved M2B and M3B phases show better linear MHD and AE stability, so this optimization trend should be considered together with an improved NBI operational regime. Indeed, similar optimization trends were already observed by other authors identifying a reduction of the energetic particle transport in configurations with low $q_{min}$, although enhanced in configuration with high $q_{min}$ [97, 98, 99, 100].

5 Destabilizing effect of the energetic particles on MHD modes

In this section we analyze the coexistence between MHD modes and low frequency AE, as well as the destabilizing effect of the energetic particles on MHD modes. If the MHD activity is dominant but the growth rate of the sub-dominant low frequency AE (BAE/BAAE) is not negligible (1/2 to 1/3 that of the dominant MHD modes), the MHD modes are further destabilized by the energetic particle drive so the MHD mode growth rate and frequency also change with the NBI injection intensity and beam energy (see Figure 14, panels a) and c), modes $n=1$ and $n=2$). There is a smooth transition in terms of growth rate and frequency between MHD modes and low frequency AEs as the NBI injection intensity increases, pointing out that both instabilities should be partially overlapped and coexist in a range of $\beta_{f}$ values. On the other hand, the growth rate and frequency of MHD/AE instabilities show a stronger dependency with the beam energy, leading to a distinct transition between MHD modes and AEs, see for example the evolution of the $n=2$ MHD mode if $V_{th,f}/V_{A0}=0.1$ to a $n=2$ BAE if $V_{th,f}/V_{A0}=0.5$ (figure 14, panels b and d).

In summary, DIII-D reverse magnetic shear discharges show an overlapping and possible coupling between MHD and AE low frequency instabilities, because the source of free energy added by the energetic particle density gradient further destabilizes unstable MHD modes in addition to the kinetic resonant destabilization of AEs. It should be noted that the destabilizing effect of the energetic particles on MHD modes is weak if the MHD modes growth rate is large compared to low frequency AE growth rate, negligible if the dominant MHD mode growth rate is 3 times larger. In addition, if the NBI operational regime is optimized to minimize AE activity (operation below the critical $\beta_{f}$ or in the weakly resonant regime), the overlapping between MHD modes and low frequency AE is weaker, optimizing the device performance with respect to low frequency AEs.

6 Conclusions and discussion

The simulation results of dominant and subdominant modes in reverse magnetic shear discharges in DIII-D are in reasonable agreement with measurements. MHD modes are unstable as soon as the reverse shear region appears in the plasma, further enhanced as the reverse shear region deepens. In addition, unstable AEs may limit the device performance at the beginning of the discharge, before the formation of the reverse shear region, remaining unstable as subdominant modes through the rest of the discharge.

MHD modes are destabilized in the bottom of the reverse shear region near the $q=2$ rational surface, although if the $q=3$ rational surface is also resonant in the inner plasma region, as in the shot $164841$, MHD modes are destabilized close to the magnetic axis. AEs are unstable near the magnetic axis in the discharge with positive/negative magnetic shear at the plasma core, shots $164842$ and $164922$, because the magnetic shear at the plasma core is weak in both discharges. Core negative shear configurations can have a larger magnetic shear compared to core positive shear configurations avoiding the presence of main rational surfaces in the plasma core and in the reverse shear region, improving MHD and AE stability.

Thanks to the analysis of the subdominant modes in the simulations we can study the stability of different AE families as overtones of the same perturbation, reproducing the wide spectrum of AE frequencies observed in the experiment. For example, the numerical model predicts unstable AEs in the frequency range of $f=[125,350]$ kHz during the phases B and C of shot $164842$, reproducing the CO2 interferometry measurements. The magnetic diagnostic shows the destabilization of low frequency AE, TAE and MHD modes that coexist in a frequency range between $f=[5,80]$ kHz. It should be noted that some discrepancies between the experimental data and the simulations are identified, for example the destabilization of modes in the frequency range of $f=70$ to $125$ kHz not observed in the experiment. Such a mismatch can be explained by the stabilizing effect of the EP finite Larmor radius or electron-ion Landau damping, not considered in the study. Also, the fact that the safety factor and energetic particle profiles are not directly measured during the experiment, and the profiles can sensitively affect the AE stability threshold.

The NBI operational regime can be optimized to reduce the AE activity if the NBI deposition is off-axis and the beam energy increases (compared to the experimental NBI operational regime), reducing the AE frequency and growth rate and stabilizing EAE/NAE. If the beam energy is further increased, the NBI operation can enter in the non resonant regime, where TAEs, low frequency AEs or RSAE are stable, as was predicted and observed in previous studies of AE stability in LHD, TJ-II and DIII-D [89, 90, 92]. It should be noted that reducing the beam energy decreases the number of lost EPs for the same injection intensity, so this optimization trend can’t be ignored. On the other hand, the NBI injection intensity is above the critical $\beta_{f}$ to destabilize AEs, so no improvement is observed in the simulations even if $\beta_{f}$ is reduced down to $0.01$. The magnetic field configuration can also be optimized because there are operation windows with improved MHD and AE linear stability, particularly if the $q$ profile in the plasma core and reverse shear region is located between the $q=1$ and $q=2$ rational surfaces.

Energetic particle drive can further destabilize MHD modes whose growth rate and frequency change with the NBI injection intensity and beam energy. The simulations show an overlapping between MHD modes and low frequency AE as well as a smooth transition between MHD modes and AE if $\beta_{f}$ increases, so cross-perturbation effects should be expected. Nevertheless, MHD modes and low frequency AE overlapping can be avoided or minimized if the NBI operates below the critical $\beta_{f}$ or in the weak resonance regime, improving DIII-D performance.

Appendix

Mode identification

Figure 15 shows the eigen-function of several instabilities obtained in the simulations. It should be noted that the AEs are primarily identified by the instability position on the Alfvén continua (see Figure 1), although some information can also be gained by analyzing the eigenfunction structure, for example from the dominant modes overlapping (overlapped consecutive poloidal modes indicates the destabilization of a TAE) or the eigenfunction peak location along the normalized minor radius (the RSAEs are destabilized near the minimum of the q profile by a single dominant mode). Figure 15a shows a $n=2$ interchange mode (MHD-like instability) destabilized near the q=3 resonant surface during the M1A phase, characterized by a dominant single mode ($-5/-2$) and an instability frequency lower than $1$ kHz. Figure 15b shows a $n=2$ BAE destabilized during the M20 phase, also dominated by a single mode although in this case both parities are important ($6/2$ and $-6/-2$) and the instability frequency is $30$ kHz. Figure 15c shows an $n=5$ TAE unstable in the M10 phase where the $13/5$ and $14/5$ modes are coupled and the instability frequency is $155$ kHz. Figure 15d shows a $n=6$ EAE destabilized during M3C phase where $10/6$ and $12/6$ modes are coupled and the instability frequency is $235$ kHz. Figure 15e shows a $n=6$ NAE destabilized during M2B phase where the modes $-15/-6$ and $-18/-6$ are coupled and the instability frequency is $500$ kHz. Figure 15f shows a $n=6$ RSAE destabilized during M1B phase, where there is a single dominant mode $11/6$ located near the minimum of the reverse shear region with a frequency of $115$ kHz.

Figure 16 shows the Electron Cyclotron emission (ECE) data and the electron temperature fluctuations ($\delta T$). The ECE data is measured along chords that span between $R=[1.4,2.2]$ m. The panels a and b show the data of the shot 164841 at $t=1650$ ms for $f=38$ kHz, indicating the presence of an instability in the range of the $30-40$ kHz located in the inner plasma region, matching the $6/2$ BAE reproduced by the simulations. The panels c and d show the data of the shot 164842 at $t=1720$ ms for $f=158$ kHz, pointing out an instability in the range of the $155-165$ kHz located in the middle-outer plasma region, similar to the $13/5-14/5$ TAE calculated by the simulations although located in the middle plasma region. The panels e and f show the data of the shot 164922 at $t=2910$ ms for $f=228$ kHz, where an instability in the range of the $225-235$ kHz is observed in the inner plasma region, congruent with the $10/6-12/6$ EAE obtained by the simulations. The NAE is not compared with the experimental data because the maximum ECE frequency is limited to $250$ kHz. For the RSAE, the experimental data shows instabilities with the characteristic frequency swiping of the RSAE although at frequencies above $140$ kHz, thus this case is neglected from the comparative study due to that discrepancy.

This material based on work is supported both by the U.S. Department of Energy, Office of Science, under Contract DE-AC05-00OR22725 with UT-Battelle, LLC and U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences, using the DIII-D National Fusion Facility, a DOE Office of Science user facility, under Award No. DE-FC02-04ER54698. This research was sponsored in part by the Ministerio of Economia y Competitividad of Spain under project no. ENE2015-68265-P. DIII-D data shown in this paper can be obtained in digital format by following the links at https://fusion.gat.com/global/D3D_DMP.

DISCLAIMER

This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.