Theory and observation of the onset of nonlinear structures due to eigenmode destabilization by fast ions in tokamaks

Supra-thermal fast ions exist in fusion-grade tokamaks as a result of neutral beam injection (NBI), resonant heating, and fusion reactions. This population of energetic particles (EPs) can strongly resonate with Alfvénic modes and excite instabilities that can seriously deteriorate the confinement.^1–5 The control of this interaction is necessary for the achievement of burning plasmas scenarios, in which the fast ions need to have sufficient time to transfer their energy to the background—mostly from drag (slowing down) on electrons—in order to ensure high temperature for the continuation of fusion reactions. This energy transfer mechanism is considered essential for the good performance of ITER.

Theoretically, the time evolution of the amplitude of a mode interacting with EPs can exhibit a variety of patterns as the mode departs from an initial linear phase towards its nonlinear response. During this evolution, several bifurcations can take place, with typical phases being steady, regular, chaotic, and chirping oscillations.⁶ Upon the kinetic interaction of particles with an eigenmode, nonlinear phase-space structures may spontaneously emerge in the resonance regions of the particle distribution, depending on the competition between drive, damping, and collisionality.⁷ These disturbances can self-consistently support anharmonic oscillations, in a generalization of Bernstein-Greene-Kruskal (BGK) modes⁸ which, in the presence of wave damping, are pushed towards lower energy states. These self-trapped structures consist of accumulation and depletion of particles in phase-space and are commonly referred to as clumps and holes, respectively.

Frequency chirping can emerge just above the threshold for marginal stability where the energy extracted from resonant EPs slightly exceeds the power being absorbed by the background dissipation, as discussed in Ref. 7. The initial nonlinear response is to relax the EP distribution in its resonance region, which would reduce the drive which can then damp the mode. However, the plasma-EP system can also find an alternate option, of slightly shifting its frequency, thereby still tapping the free energy of previously untapped neighboring non-resonant particles that then become resonant due to the changed frequency. In the fully developed chirping state, the resonant regions of an enhanced number of particles (clumps) or of a deficient number of particles (holes) are trapped by the finite amplitude wave, and these regions of phase space shift the resonant distribution to lower energy regions of phase space, with the released energy being absorbed by the background dissipation mechanisms that are present while keeping the nonlinear amplitude hardly changed. Therefore, the frequency variation itself allows for the moving structures to access phase space regions with distribution function gradients otherwise inaccessible, which leads to convective losses over an extended region.

For the sake of clarity, we distinguish between the terminologies, frequency chirping and frequency sweeping, which often appear in the literature. While the former is normally associated with the fast kinetic response of resonant particles (typically of order 1 ms) studied in this work, the latter usually corresponds to a slow evolution (typically 100 ms) of Alfvénic modes, in the presence of a non-monotonic safety factor q profile that relies on the time variation of q_min. Sweeping events are associated with a modification of the plasma equilibrium (and consequently, of the dispersion relation), for example, in the case of Alfvén Cascades.⁹ Chirping is faster and harder to suppress using external control.

Chirping modes can have frequency shifts greater than the linear growth rate γ_L and are observed to be precursors of even worse scenarios, known as avalanches. A spectrogram showing repetitive chirping cycles followed by avalanches for toroidal Alfvén eigenmodes (TAEs) in National Spherical Torus Experiment (NSTX), for several toroidal mode numbers, is presented in Fig. 1(a). The inset shows four of the chirping events and indicates that it consists mostly of a down-chirping. The system preference for a direction (up or down in frequency) has been theoretically linked with the competition between different collisional processes.¹⁰ Figure 1(b) shows very significant neutron rate drop correlating with the avalanches.

The long-range chirping evolution was described by the Berk-Breizman prediction for the frequency variation $δω$ scaling with the nonlinear bounce frequency ω_b to the power of 3/2.⁷ It has been successfully used for applications that include the inference of mode amplitude on MAST¹¹ and the estimation of kinetic parameters such as drive and damping in JT-60U¹² and in NSTX.¹³ 

The present work however focuses on establishing the conditions for chirping onset rather than modeling their long-term evolution, in order to predict the likely character of EP transport. The EP losses are typically diffusive (e.g., due to mode overlap, RF fields, turbulence, and collisional scattering) or convective (as a result of chirping oscillations and collisional drag). We describe the methodology for the generalization of previous works and include micro-turbulent stochasticity, which [for the case of no ion cyclotron resonant heating (ICRH)] has been shown to compete with, and even greatly exceed, collisional scattering in many tokamak experiments¹⁴ and therefore needs to be added to the stochasticity due to pitch-angle scattering. We note that RF-induced diffusion of resonant ions can play a similar role in altering phase-space chirping structures.^15–18 In this work, however, we do not examine experimental cases in which the RF decorrelation mechanism is dominant.

Chirping and quasilinear (QL) regimes correspond to two distinct limits of kinetic theory. Since they may be competing mechanisms in the modification of the distribution of fast ions in tokamaks (and their consequent transport), their parameter-space regions of applicability need to be carefully addressed. The derivation of the QL diffusion equations^19,20 relies on averages, over a statistical ensemble, which smooth out the distribution function. In order to justify the resulting smooth, coarse-grained distribution, stochastic processes (which can be intrinsic due to mode overlap or extrinsic due to collisions) need to be invoked. The fast-varying response associated with the ballistic term is disregarded, which implies that entropy is no longer conserved. Consequently, QL theory kills phase correlations and cannot capture chirping events, since chirping needs time coherence from one bounce to the next in order to move nonlinear structures altogether over phase space. QL diffusion theory needs phase decorrelation, i.e., particles need to be expelled from a phase-space resonant island at a time less than the nonlinear bounce time. This means that there are no particles effectively trapped. Due to the reduced dimensionality of phase space, the QL description is less computationally demanding than the full nonlinear description needed to capture chirping. It is also much less computationally expensive than particle codes. A criterion for chirping likelihood is an important element for identification of parameter space for QL theory applicability for practical cases and consequent validation of reduced models. An example of such models is the Resonance-Broadened Quasilinear (RBQ) code.^21,22 It uses the usual structure of the QL system written in action-angle variables²³ with a broadened resonance width that scales with bounce frequency, growth rate, and collisional frequency.^24,25 In this work, we build predictive capabilities regarding the likelihood of the nonlinear regime, which can be useful for burning plasma scenarios. If further validated and verified, the developed methodologies could be of practical importance for predictive tools of EP distribution relaxations in the presence of Alfvénic instabilities.

This paper builds on ideas introduced in Ref. 26 and presents new and detailed comparison of chirping likelihood for modes measured in NSTX, DIII-D, and Tokamak Fusion Test Reactor (TFTR). The detailed analysis was made possible by the numerical implementation of novel tools into the framework (e.g., for the phase-space resonance averaging procedure and for the mode structure experimental comparison), which are described here. The theoretical framework is extended and properties of the governing equations for the chirping prediction are unveiled. In addition, an original study of how the value of the particle resonant speed itself affects the wave oscillatory nature is presented, which is shown to be in agreement with observations on JT-60U. This paper is organized as follows. In Sec. II, the proposed theoretical methodology is presented. In Sec. III, the numerical procedure is presented along with analyses of experimental results. Discussions are presented in Sec. IV, and the  Appendix is devoted to discussions on the chirping likelihood in terms of the beam injection energy.

In axisymmetric tokamaks, $E, Pφ$⁠, and μ (corresponding to energy, canonical toroidal angular momentum, and magnetic moment, respectively) are considered invariants of the unperturbed motion for EPs interacting with modes that have frequencies much lower than the cyclotron frequency. Their expressions, all per unit mass of EPs, are given in S.I. units by

where m_EP and q_EP are the mass and charge of EPs, v is the EP speed, R is the tokamak major radius, ψ is the poloidal flux divided by $2π$⁠, B is the magnitude of the magnetic field, and θ and $φ$ are the poloidal and toroidal angles.

If the wave-particle interaction Hamiltonian is assumed, like in the NOVA code,²⁷ to have dependences on time and poloidal and toroidal angles as $ei(mθ−nφ−ωt)$⁠, where m and n are the wave numbers associated with them, it follows that upon particle interaction with a mode, the invariances of $E$ and $Pφ$ are broken but a new invariant arises: $E′=E+(ω/n)Pφ$ while μ is kept nearly constant. The existence of two invariants implies that the resonant particle dynamics is essentially one-dimensional.²⁸ Therefore, a new variable $I=−Pφ/n$ (at constant $E′$⁠) can be used to describe this relevant one-dimensional path for EP dynamics (for steepest distribution function f modification). The projection of the gradient operator onto this path is then given by

The resonances for the several harmonics, which are multiple surfaces in $(E,Pφ,μ)$ space, are defined by

where j is an integer and $ωφ$ and $ωθ$ are the local toroidal and poloidal transit frequencies, respectively. The phase-space integration is represented as

where c is the light speed and $σ∥$ accounts for counter- and co-passing particles.

Collisional processes are an important element in the determination of the nonlinear character of wave oscillations.^6,28,29 Stochastic processes, such as pitch-angle scattering, act to destroy the coherent structures that support wave chirping while drag, or slowing down, is formally equivalent to chirping and enhances the convective transport of these nonlinear structures. For EPs, the Fokker-Planck collisional operator that enters the kinetic equation can be approximately written as a superposition of pitch-angle scattering and drag as follows:³⁰ 

where f is the distribution function; $ν⊥$ and $τs−1$ are the $90$° pitch angle scattering rate and the inverse slowing down time, which are given by

The critical speed, above which electron drag dominates over ion drag, is defined as

and

For a given ion species β, $Aβ=mβ/mp$⁠, with m_p being the proton mass, $lnΛ$ is the Coulomb logarithm, Z_i and Z_EP are the background and energetic ion atomic number, n_a is the density, T_a is the temperature, m_a is the mass, and $vTa=2Ta/ma$ is the thermal speed for a given species a. The three-dimensional velocity derivatives are projected onto the aforementioned path of the steepest gradient of the distribution function so that

Since $∂Pφ∂v=Rφ̂$ and $vφ≈v∥BφB$⁠, it follows that the effective collisional operator for fast ions resonating with a given mode can be approximately cast in the form

where the effective scattering ν_scatt and drag ν_drag coefficients are

and

In order to be properly evaluated, the above expressions need to be averaged over the orbit bounce motion of a particle, as detailed in Sec. III C.

Stochastic diffusion is determined by collisional scattering processes, such as pitch angle scattering, as well as additional processes, such as the effect of the background turbulence that can be dominant in the determination of the global heat outflow. Micro-turbulence is introduced using a procedure that follows the one introduced in the pioneering work by Lang and Fu,¹⁴ where it was considered that the main contribution comes from electrostatic ion temperature gradient (ITG) turbulence via radial diffusion rather than the velocity diffusion. As with collisional scattering, the turbulent diffusion operator is projected onto the relevant one-dimensional path for particle dynamics, represented by the variable Ω. In the limit of a large aspect ratio, circular cross-section torus, it is possible to make use of the relation $∂∂r=∂Ω∂Pφ|E′∂Pφ∂r∂∂Ω$ with $∂Pφ∂r≈−qEPmEP∂ψ∂r$⁠, and the spatial micro-turbulence diffusion operator along the radial coordinate r can be written as

where D_EP is the EP diffusivity. Therefore, the ratio between the effective stochasticities coming from micro-turbulent [Eq. (8)] and collisional processes [Eq. (6)] is

A similar ratio was derived in Ref. 14, where a Lorentz operator was used for the pitch-angle scattering.

A subtlety in applying the model relies on the determination of D_EP. It is expected to be lower than thermal ion diffusivity $Dth,i$ since the micro-turbulence wavelength is typically smaller than the beam cyclotron orbit. Historically, the effect of micro-turbulence on EP transport has been neglected based on the fact that since EPs have large orbits, they should experience several phases of the turbulent fields in such a way as to cancel out its overall effect. Although fast ion turbulent transport is negligible compared to Alfvénic-induced transport (see, for example, studies on ASDEX-U³¹ and DIII-D³²), turbulence can be an important transport mechanism, as compared to collisions, at the onset of the evolution of modes, when their amplitude is still small. Over the past decade, the modeling of D_EP has been studied by several groups.^33–41 Following Ref. 14, we have chosen to determine D_EP from the scalings that follow from gyrokinetic simulations of Gyrokinetic Toroidal Code (GTC),³⁸ in which D_EP is proportional to the thermal ion diffusivity, D_i, and a function of $Ti/EEP$⁠, where E_EP is the energy of the EPs. The GTC simulations assumed a specific plasma background that can be considerably different from a given discharge being analyzed. Therefore, a significant error can be expected to be associated with the inferred value of D_EP. Alternative ways of obtaining D_EP could be achieved by using Refs. 33 and 34. In our analysis, we infer the value of D_i from the outputs of the global transport code TRANSP^30,42 at the position where the mode structure is peaked at the time being analyzed. The particle diffusivity is known to have a huge error associated with it because TRANSP cannot resolve well particle sources, especially close to the wall. On the other hand, the thermal conductivity χ is reasonably well known and therefore is frequently used⁴³ as an indication of the actual value of D (the exact relation for a Maxwellian distribution would be $D=2χ/3$⁠).

The onset of a mode amplitude evolution can be studied using perturbation theory (assuming small deviations of particles from their unperturbed orbits) within the kinetic framework considering that the system is close to marginal stability. References 6, 28, 29, and 44 showed that, for $ωb≪|γL−γd|$ (growth rate minus damping rate), truncation of mode amplitude at third order is justifiable. Taking ν_stoch (the overall stochasticity felt by EPs, which includes ν_scatt) and ν_drag independent of time but dependent on phase space localization, the equation for the early time perturbative mode (a mode that exists without accounting for the kinetic component) amplitude evolution can be written as a time-delayed, integro-differential cubic equation

where $H=2πωδ(Ωj)|Vn,j|4(∂Ωj∂I)3∂f∂Ω$ and

or

accounts for the wave-particle energy exchange, where e is the electric field eigenstructure and v is the velocity of a resonant particle. $l1,l2,l3$ are the integers and $ξ1,ξ2,ξ3$ are the angles conjugated to the invariants of motion (actions of the Hamiltonian). In Eq. (10), the circumflex denotes normalization with respect to $γ=γL−γd$ and t is the time normalized with the same quantity. A is the normalized complex mode amplitude of an eigenmode oscillating with frequency ω.

The phase factor $ϕ$ is a measure of the non-perturbative aspect of the linear problem. Here, we limit ourselves to the case $ϕ=2πh,;h∈ℤ$⁠, which is the case where a positive energy mode exists in the absence of the kinetic response of EPs. Then, the kinetic response leads to a perturbative response of the positive energy wave to the kinetic interaction. Note that the case $ϕ=π+2πh$ leads to a perturbative response of a negative energy wave to the kinetic interaction, while other values of $ϕ$ imply a non-perturbative response, where there is no linear wave present in the absence of the EPs. In other words, $ϕ$ represents the ratio between the dissipative and the reactive responses of the perturbing field.

The solutions of Eq. (10) can exhibit several bifurcations and therefore several phases, as shown for a bump-on-tail configuration in Ref. 6. Interestingly, Eq. (10) allows for the excitation of sub-critical instabilities and for nonlinear frequency splitting.¹⁸ If the nonlinearity in Eq. (10) is weak, the system will most likely saturate close to the linear stability state, where the trapping frequency ω_b satisfies $ωb≪γL≃γd$⁠. However, in case the solution of the cubic equation explodes, the system enters a strong nonlinear phase, which may lead to chirping modes. Indeed, long-range numerical simulations indicate the explosive behavior of A as the precursor to the formation of hole and clump structures.⁴⁵ Furthermore, chirping events are significantly enhanced by the coherence introduced by dynamical friction (i.e.,, particle drag)^10,29 and are inhibited by stochasticity from diffusive processes,¹⁵ such as resonant particle heating and collisional pitch-angle scattering, and from background turbulence, all of which contribute to causing particles to detune from a resonance. Stochastic events lead to the loss of phase information which contribute to destroy coherent structures.

Equation (10) was originally derived for a bump-on-tail system with Krook collisions⁶ and later generalized to complex tokamak geometries and also to include collisional scattering.²⁸ Lilley et al.²⁹ included the effects of drag on the bump-on-tail cubic equation and derived a criterion to determine stable and unstable regions of solutions of the cubic equation in drag vs scattering collisional parameter space. Our work aims at improving this prediction by using realistic resonant surfaces and mode structure information, coming from the NOVA and NOVA-K codes. To this end, orbit and phase space averages are employed in order to account for the effective Fokker-Planck collisional coefficients. Experimental data are analyzed in order to verify whether chirping events coincide with the occurrence with the “explosive” phase of the cubic equation, as predicted by the theory.

It has been recently shown²⁶ that a simplified bump-on-tail approach that only accounts for a single representative value for the collisional coefficients is insufficient to make predictions for a mode nonlinear nature in tokamaks. The missing physics were shown to be the absence of non-uniform mode structures, (multiple) resonance surfaces, and poloidal bounce averages that account for particle trajectories on a poloidal cross section.

A necessary but not sufficient condition for the existence of fixed-frequency, steady-state solutions in the form $A(t)=|A0|eibt$⁠, with A₀ and b constants, is that the real component of the right-hand side of Eq. (10) should be negative at late times (⁠$t→∞$⁠) when the response is stationary, i.e., when the nonlinear term is allowed to balance the linear growth

Note that the dependence of the arguments of the exponentials on b cancels out in all three terms of Eq. (10). Next, the integration over τ₁ can be straighforwardly performed in Eq. (12). Furthermore, the delta function $δ(Ωj(Pφ,E,μ))$⁠, contained in the expression for $H$⁠, restricts the integration over phase space [as defined by Eq. (2)], allowing the resonance condition to be exploited to eliminate the integral over energy space. Then, upon a re-normalization of the time variable to $z=ν̂dragτ$⁠, the following criterion for the existence of fixed-frequency oscillations is obtained²⁶ 

where

and N is a normalization for Crt, which consists in the same sums and integrations that appear in the numerator of (13) but in the absence of Int. In Eqs. (13) and (14), the quantities τ_b, ν_drag, ν_stoch, $Vn,j$⁠, and Ω_j are understood to be evaluated at $E=E′−ωPφ/n$⁠. Crt involves a phase-space integral restricted to the resonance surfaces. Although Crt does not explicitly depend on the growth rate, it does depend on the local phase-space gradient of the distribution function through the term $∂f/∂I$⁠. Because of the Crt dependence on $Vn,j$⁠, the regions where the mode eigenstructure peaks contribute most to Crt, more specifically, in regions where the resonant particle current maximizes its projection onto the eigenmode electric field. Criterion (13) was incorporated into NOVA-K making use of a polynomial interpolation for Int. Crt provides a prediction for the likelihood of a fully nonlinear phenomenon obtained only from pure linear physics elements and therefore can be tested in linear codes. This is a considerable advantage in efficiency for making a prediction of a nonlinear property. The integrand of Int is plotted in Fig. (2) for different values of $νstoch/νdrag$⁠.

Drag introduces an oscillatory behavior of the integrand of Int. This has an effect of flipping the sign of the integral kernel of the cubic equation. This causes dA/dt to vary more abruptly and therefore prevents a steady state solution from being achieved. Therefore, nonlinear chirping solution is more likely to be achieved in drag-dominated regimes. It is interesting to note that formally, drag enters the kinetic equation in a mathematically similar way as a chirping frequency does, with $kvνdrag$ replaced by $dω/dt$⁠.

Int, which is a function of phase space, is plotted in Fig. 3 as a function of $νstoch/νdrag$⁠. The positive and negative domains of Int are roughly an order of magnitude different in the interval $0≤νstoch/νdrag≲2$⁠, while for $νstoch/νdrag≫1$⁠, $Int≈1.022(νstoch/νdrag)−4$⁠.

For our investigation, we focus on the criterion for the existence or non-existence of a fixed-frequency, steady solution. To determine the physical parameters needed in the calculation described below, the global transport code TRANSP^30,42 is used. The code determines the needed particle diffusivity that matches the observed plasma parameters to the particle and heat sources. We interpret this empirically obtained result as due to both classical (which includes neoclassical) transport and micro-turbulence. The resulting thermal particle diffusivity is then used together with a model to scale the EP diffusivity given the turbulent background diffusivity, to estimate the diffusivity of EPs, which will then contribute to the value of ν_stoch. Here, we perform a quantitative study based on the time delayed cubic equation using experimental results to determine whether the theoretically predicted nonlinear character of the response correlates with the observation.

NOVA²⁷ is a nonvariational ideal MHD code primarily used to integrate non-Hermitian integro-differential eigenmode equations in the presence of EPs, using a general flux coordinate system. It uses realistic numerically calculated MHD equilibria and hence it is suited to study both conventional and spherical tokamaks. The code uses Fourier expansion in θ and cubic spline finite elements in the radial ψ direction. NOVA provides the eigenstructures used in our analysis.

NOVA-K^46,47 takes into account finite Larmor radius and orbit width effects to study the destabilization of MHD modes from the EPs free energy. We use NOVA-K to calculate the resonance surfaces, in $(Pφ,E,μ)$ space, associated with the modes. These surfaces are needed to calculate the growth and damping rates of eigenmodes in the presence of EPs as well as the integral Int. TRANSP is used to provide the distribution function which contains the necessary information on the most representative EP population as input for NOVA-K runs.

The quadratic form of MHD is particularly useful when stability analysis is addressed. If the mode frequency is assumed to be $ω+iγ$⁠, with $|γ|≪|ω|$⁠, where ω is the mode eigenfrequency and γ is its growth rate, we obtain^27,46

where $δK=ω2∫ρ|ξ|2dr/2$ is the inertial (kinetic) energy, $δWk$ is the potential energy associated with the nonadiabatic component of the distribution function, and $ωδWk$ is the power released by the resonant particles. The growth rate can be expressed as

with the diamagnetic frequency being defined by

G_mj represents mode structure matrix elements (as defined in Sec. II D of Ref. 27), which are associated with the projection of the resonant particle current onto the wave electric field. For each isolated resonance, the reduced Hamiltonian can be written as²⁸ 

where $ζ≡l1ξ1+l2ξ2+l3ξ3$ and the bounce frequency of the most deeply trapped particles can be shown to be

where the subscript r denotes the resonance location and the derivative is defined by Eq. (1). The G matrices are related to V [given by Eq. (11)] via

In order to characterize the mode being observed in the experiment, NSTX and TFTR reflectometer measurements are compared to the mode structures computed by NOVA. Reflectometry diagnostic measures the density fluctuation of the plasma at the location where the launched wave has a cutoff. The fluid displacement times the local density gradient is equivalent to the density fluctuation. All poloidal harmonics calculated by NOVA are summed up for this analysis. Modes are categorized according to their mode structure and whether their frequencies fall into a given gap of the continuum. The reflectometer cannot resolve the density in the core for the cases of the flat density profile and especially when the peak of the density is displaced from the plasma center, as typically observed during the H-mode regime. In DIII-D, Electron Cyclotron Emission is used for the purpose of measuring the mode structure, following a methodology described in Ref. 48. An interesting aspect of the observation is that, in spite of the intrinsic nonlinear nature of chirping events, the structure of the chirping mode is not substantially changed during the time evolution of the system, in accordance with what was previously reported in Refs. 13, 49, and 50. This gives credence to the use of a linear code for the eigenstructure identification at early times.

In order to calculate the above expressions for the effective collisional coefficients (6) and (7), NOVA-K is employed to perform bounce averaged calculations since the bounce motion is much faster than the perturbative mode evolution. Then, a phase-space average needs to be performed to account for the contribution of resonance surfaces spread over phase space, as described below.

The period of a particle poloidal bounce or transit motion is numerically calculated following the framework of Ref. 27. The formulation captures both finite orbit width (FOW) and finite Larmor radius (FLR) effects and accounts for particle drifts due to curvature and $∇B$ terms. The averaging procedure is taken via Bessel functions. For each set $(E,Pφ,μ)$⁠, the trajectories are a priori known, and hence, there is no additional need to explicitly follow the particle motion within the code.

The average over phase space is taken along the surfaces over which the resonance condition is satisfied, for different poloidal bounce harmonics. The phase space volume elements are weighted in accord with their relative contribution to the overall growth rate. Specifically, we evaluate

where $Q=1ω∑j〈qEPe·v〉2∂f∂Iδ(Ωj)$ is the contribution to the growth rate, γ_L, from a given phase space location that satisfies $Ωj(E,Pφ,μ)=0$⁠. In the present study, the phase-space averages are taken over several harmonics of a given mode. This averaging technique was previously used to predict the TAE amplitude saturation in TFTR experiments.⁴⁷ 

In order to show the importance of micro-turbulence as a chirping suppression mechanism, an n = 4 TAE driven by alpha particles in the Tokamak Fusion Test Reactor (TFTR) shot 103101 was studied. This mode, which was observed to oscillate at a constant frequency, was analyzed in terms of the proposed criterion discussed above. Its frequency and strength relative to the background field are presented in Fig. 4. The corresponding mode structure obtained with the NOVA code is shown in Fig. 5. Shown in Fig. 6 are scans over the inferred experimental values of ν_stoch and ν_drag, denoted by $νstoch(exp)$ and $νdrag(exp)$⁠, for the situations (a) without the inclusion of micro-turbulent stochasticity and (b) with its inclusion.

A detailed visualization of how far the experiment is from the expected boundary that separates the steady and chirping regions, as predicted by the model, can be provided by scanning the criterion integral for several values multiplying the actual ν_drag and ν_stoch (as shown in Fig. 6). The vertical axis represents the constants that are multiplying ν_stoch, while the horizontal one represents constants that are multiplying ν_drag. The point $(νstoch/νstoch(exp)=1,νdrag/νdrag(exp)=1)$ corresponds to the inferred experimental condition, and the criterion boundary is represented by the 0.0 contour. Contours of negative values represent regions in which chirping is expected, while the positive ones represent expected steady-state (fixed frequency). If experimental conditions imply a positive Crt sufficiently far from this transition boundary, steady solutions should arise in experiment.

Each contour plot corresponds to a given value of Crt [Eq. (13)], as labeled in the two plots in Fig. 6. In part (a), where Crt is evaluated only in terms of collisional pitch-angle scattering and drag, Crt is a negative number. Upon the addition of micro-turbulent scattering, the re-evaluation of Crt shows that the point (1, 1) indicating the experimental conditions now operates where Crt > 0, which means that the criterion predicts that the mode should not be able to chirp, being in agreement with the observation. For this case, the estimation for D_EP based on the scaling of Ref. 38 led to $DEP≈0.1 m2/s$⁠, which is consistent with direct measurements using beam blips in similar experiments carried out in TFTR.⁵¹ Numerically calculated $νscatt(exp)$ and $νdrag(exp)$ were $8568 s−1$ and $3138 s−1$⁠, respectively. $νstoch(exp)$ was found to be nearly 10 times larger than $νscatt(exp)$ for this particular mode. The $90$° pitch-angle scattering and the inverse slowing down time were found to be $ν⊥=0.19 s−1$ and $τs−1=4.55 s−1$⁠.

The NOVA-K code does not account for turbulence stochasticity, and consequently, it evaluates the criterion via an indirect method, using the following procedure. First, we note that Crt only depends on the ratio of the overall stochasticity and collisional drag. From Eqs. (3), (4), (6), and (7), we see that the effective drag coefficient is proportional to $Te−3/2$ while the effective collisional scattering is insensitive to T_e. It is then possible to mimic the inclusion of turbulence into ν_stoch by artificially decreasing drag, which can be achieved by multiplying T_e by a corresponding numerical factor that leads to the correct value of the averaged values of $νstoch/νdrag$⁠.

A series of dedicated shots were performed on DIII-D with the objective of triggering chirping and studying what conditions most strongly determine a mode nonlinear evolution into the chirping regime. These shots used high ion temperature at the core (⁠$10−12 keV$⁠), $qmin ∼2$⁠, and strong toroidal rotation (up to 50 kHz on axis). A good example is the chirping observed in shot 152828, shown in Figs. 7 and 8.

The frequency of the chirping mode was too low to be a TAE. At first, NOVA was run in a mode that enabled only the Alfvénic branch to be captured and no reasonable mode structure was found. Then, NOVA was run allowing both the Alfvénic and the sonic branches. A mode that matches experimental evidence was obtained and identified to have the characteristics of a beta-induced Alfvén-acoustic eigenmode (BAAE)⁵³ and is shown in Fig. 9. The mode was identified using the data available on mode localization from the electron cyclotron emission (ECE), frequency, local rotation, and toroidal mode number (n = 1 for the chirping mode). For the comparison with ECE data for electron temperature fluctuation, the fluid displacement $ξ$ calculated by NOVA was post-processed using the relation

where in this equation $γ=5/3$ is the ratio between specific heats.

The time-delayed cubic Eq. (10) is derived assuming that the nonlinear bounce frequency ω_b is much less than the net growth rate. This implies a limitation on the cubic equation, with its domain of validity being restricted to the early nonlinear phase, when mode amplitude, represented by ω_b, is still small compared to γ. This framework therefore cannot be used to model the full chirping events. From DIII-D observations, we have noted that the chirping behavior typically starts when χ_i drops to values lower than $0.25 m2/s$²⁶ but chirping persists even when χ_i raises to values that would not admit the onset of the chirping process, according to the theory. The chirping cycle appears to involve some degree of hysteresis since once the first chirp happens, there is a tendency of the system to continue chirping. Thus, this observation appears to indicate that once the chirping structures are already embedded in the system, continued new chirping can still arise even though the apparent diffusivity has increased to the point that chirping would not occur if phase space structures were not already present. This hysteresis can be noted in Fig. 8 as the chirping persists up to the point where $χi≈1.5 m2/s$⁠.

Chirping does not restart when χ_i drops for the second time, at $t=1.05 s$⁠. This is probably because of the effect of strong neoclassical tearing modes (NTM) that are present at this time (as can be seen from the parallel horizontal spectral lines that begin at around $t=992 ms$ in Fig. 8). The suppression then probably occurs because the NTMs are detuning the EPs from the resonance and therefore contributes to extinguishing the drive (Alfvén waves are not observed in Fig. 8 after $t=1 s$⁠).

Figure 10(a) shows the experimental condition for the mode in DIII-D shot 152828 shown in Fig. 9 before chirping starts (at $t=900 ms$⁠, when $Dth,i≈0.9 m2/s$⁠) and Fig. 10(b) shows it during chirping (at $t=960 ms$⁠, when $Dth,i≈0.25 m2/s$⁠), when the point has transitioned from the positive to the negative region, in agreement with the generalized criterion prediction. The criterion evaluated close to the onset of chirping, at $t=960 ms$⁠, used the mode structure calculated at $t=970 ms$ (Fig. 9) since trustworthy equilibria and profiles could not be made exactly at the time of its onset.

Similar framework of analysis was also applied to NSTX shot 141711, whose spectrogram (Fig. 1) showed chirping for several TAEs with different toroidal mode numbers. A number of chirping modes were studied, with the evaluation of Crt for all of them leading to negative numbers. A representative example is a TAE that was catalogued by comparing the mode structures calculated by NOVA with the reflectometer measurements, as shown in Fig. 11. The contour plot for Crt for this mode is presented in Fig. 12. During the chirping time window in NSTX shot 141711, the ion transport was dominated by neoclassical processes. It was observed that micro-turbulence was not strong enough to suppress the phase-space structures that sustain chirping since the (1, 1) barely moves in Fig. 12 upon the addition of turbulent stochasticity.

We have performed a study of the early phase of chirping events in tokamak plasmas by means of realistic calculations of eigenstructures and collisional coefficients. The proposed methodology has been described in fair detail. In order to generalize the theoretical predictions from the previous bump-on-tail approach, we have employed an action-angle formalism, with a similar perturbative procedure used in Ref. 28. It leads to the generalized criterion for the likelihood of chirping and fixed-frequency solutions. It should be noted that it only allows a steady solution to exist but does not guarantee that the system will necessarily evolve to such a state. It may happen that the evolution to steady states is only accessible when higher-order nonlinear terms are taken into account, which are not captured by the lowest nonlinear order perturbative approach used by the framework of the time-delayed cubic equation [Eq. (10)].

An interesting effect unveiled by the methodology used in this paper is that the disparity of magnitudes of the positive and negative regions of Int [Eq. (14)] typically leads to the prediction that a given mode should chirp, even for $〈νstoch/νdrag〉>1$⁠, unless a strong diffusive mechanism implies $〈νstoch/νdrag〉≫1$⁠. This additional mechanism is observed to be often related to turbulent processes. The contribution from regions of low μ turns out to be quite important since the magnitude of Int is larger for $νstoch/νdrag≪1$ (with ν_scatt being proportional to μ). This shows the importance of using a weighted average of the transport coefficients and not only a single representative point.

The employed approach is applicable for the modes near the stability threshold. We implement it perturbatively using the mode structures computed by the ideal MHD code NOVA. The perturbative approach is justified when the growth rates evaluated for TAEs and RSAEs in NSTX, DIII-D, and TFTR plasmas are typically much smaller than the eigenfrequencies. The case of unstable BAAEs requires more careful stability analysis that is beyond the scope of this paper. Nevertheless observations of those modes in DIII-D shown in Fig. 8(b) indicate that these modes are slowly growing, suggesting that they are near the instability threshold.

Background turbulence was introduced to the model through the estimates of EP diffusivity, which adds to the diffusivity associated with collisional scattering and contributes to the suppression of chirping. We have set up plausible rules for determining the effective diffusivity and drag. The main uncertainty in the diffusivity chosen for the EPs is the value of the diffusivity contribution from the background fluctuations. We have chosen to base this value to be in accordance with what TRANSP predicts to be the diffusion needed to produce the observed energy confinement time. However, an additional uncertainty still remains in determining how this thermal diffusion coefficient would be related to the effective diffusion coefficient for the large orbit EP. We chose to use an empirical formula obtained in Ref. 38, based on electrostatic turbulence. However, the value chosen is rough, and the actual scaling for a given experiment might differ considerably from our choice.

A number of factors that may influence the chirping formation are not captured by the current theory. For example, static 3D fields⁵⁴ have been shown to affect bursting Alfvén modes and reduce chirping. Besides that, effects such as toroidal field ripples, energy diffusion, radio frequency heating fields,^15,17 electromagnetic turbulence, and neoclassical tearing modes can be important in some scenarios to prevent chirping formation due to randomization of phase space, with consequent resonance detuning. Other limitations are that the cubic equation assumes small mode amplitude, which is not necessarily the case in the experiment, and also that no mode overlap has taken place.

Chirping events require self-trapped resonant particles to remain locked with the excited chirping frequency for successive wave-trapping bounce times. The maintenance of this time-dependent resonance condition can induce significant convective transport over an extended region of phase space. The present work can be helpful in addressing the likelihood of convective and diffusive transport of fast ions and therefore be useful as a predictive tool for present-day and next-generation devices.

We acknowledge fruitful discussions with G.-Y. Fu, E. D. Fredrickson, and A. Bierwage. This work was supported by the São Paulo Research Foundation (FAPESP, Brazil) under Grant Nos. 2012/22830-2 and 2014/03289-4 and by U.S. Department of Energy (DOE) under Contract Nos. DE-AC02-09CH11466 and DE-FC02-04ER54698. This work was carried out under the auspices of the University of São Paulo—Princeton University Partnership, project “Unveiling Efficient Ways to Relax Energetic Particle Profiles due to Alfvénic Eigenmodes in Burning Plasmas.”

In planning and interpreting experiments, it can be meaningful to understand the likelihood of observing wave chirping in terms of the resonant particle speed. This is because the injected NB ions can be either supra- or sub-Alfvénic which, for the case of TAEs, would lead to dominant resonances located around $v∥=vA$ and $v∥=vA/3$⁠, respectively. In particular, a puzzling observation in JT-60U is that the abrupt large events (ALEs) and their associated bursts and chirps usually happen when the beam is injected supra-Alfvénically, using negative-ion-based neutral beams.⁵⁵ The ALEs have been intensively studied over the past decade.^56–61 Here, we propose an explanation for the likelihood of these events in terms of the framework adopted in this paper and we show that it qualitatively agrees with JT-60U observations.

Let us analyze the case of a TAE for which the resonance condition is $Ω=ωTAE+nv∥/R−jv∥/qR=0$⁠. For an initial analysis, let us consider scenarios in which micro-turbulence can be ignored in comparison with collisional scattering. Then, the relevant ratio for the criterion, Eq. (13), would be simply $(νscatt/νdrag)3$⁠, with the effective scattering given by Eq. (6)

where ω_c is the EP cyclotron frequency and $dPφ≃−(ωc/q)rdr$⁠. The effective drag is given by (7)

with v_c given by (5). Therefore

where

and θ_vB is the angle between v and B. The ratio between the chirping-criterion-relevant parameters $(νscatt/νdrag)3$ prescribed by Eq. (A1) for the cases in which the main resonance of a TAE is at $v∥=vA$ and $v∥=vA/3$ is plotted in Fig. 13. The ratio $(νscatt/νdrag)v∥=vA3/(νscatt/νdrag)v∥=vA/33$ is 1/9 for $vc/vA≈0$⁠. If $vc cos θvB/vA≲0.5$⁠, it means that supra-Alfvénic injection implies more likelihood for chirping. It is a direct implication of the chirping criterion [Eq. (13)] and the form of Int (Fig. 3). Considering deuterium as the only background ion species, $ne≈ni$⁠, $lnΛi≈lnΛe$⁠, and using typical parameters for JT-60U, such as $B=3T, ni=3×1019m−3$⁠, and $Te=3 keV$⁠, one obtains $vcvA≈0.3$⁠. The factor $cos θvB$ contributes to reduce the variable $vc cos θvB/vA$ used in Fig. 13. Therefore, because of the resonant velocity dependence, if all other parameters are held fixed, supra-Alfvénic injection implies more chances of observing wave chirping for TAEs in JT-60U than in the sub-Alfvénic situation, provided that the mode is the same in the two situations. A similar conclusion would also hold for other conventional tokamaks, such as DIII-D.

The previous conclusion was based on the fact that the stochasticity coming from micro-turbulence is low compared to the one coming from pitch-angle scattering. However, the inclusion of micro-turbulence using the scalings found in Ref. 38 does not change the conclusion that higher resonant velocities imply more likelihood for chirping, provided that $vc/vA$ is sufficiently small. In fact, turbulence even contributes to strengthen this interpretation. This is because³⁸ it is found that the micro-turbulence diffusivity goes as $v−2$ for passing particles and $v−4$ for trapped particles, which is, similar to the effective pitch-angle scattering, an inverse velocity dependence. Therefore, micro-turbulence contributes to further reduce the ratio $(νscatt/νdrag)v∥=vA3/(νscatt/νdrag)v∥=vA/33$⁠, which makes chirping for supra-Alfvénic injection even more probable than in the purely collisional case.