# Onset of nonlinear structures due to eigenmode destabilization in tokamak plasmas

###### Abstract

A general methodology is proposed to differentiate the likelihood of energetic-particle-driven instabilities to produce frequency chirping or fixed-frequency oscillations. The method employs numerically calculated eigenstructures and multiple resonance surfaces of a given mode in the presence of energetic ion drag and stochasticity (due to collisions and micro-turbulence). Toroidicity-induced, reversed-shear and beta-induced Alfvén-acoustic eigenmodes are used as examples. Waves measured in experiments are characterized and compatibility is found between the proposed criterion predictions and the experimental observation or lack of observation of chirping behavior of Alfvénic modes in different tokamaks. It is found that the stochastic diffusion due to micro-turbulence can be the dominant energetic particle detuning mechanism near the resonances in many plasma experiments, and its strength is the key as to whether chirping solutions are likely to arise. The proposed criterion constitutes a useful predictive tool in assessing whether the nature of the transport for fast ion losses in fusion devices will be dominated by convective or diffusive processes.

## I Introduction

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 damage the confinement Heidbrink (2008); Gorelenkov et al. (2014); Lauber (2013); Chen and Zonca (2016). 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 take place, with typical phases being steady, regular, chaotic and chirping oscillations Berk et al. (1996). 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 Berk et al. (1997a). These disturbances can self-consistently support anharmonic oscillations, in a generalization of BGK modes Bernstein et al. (1957) 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. Berk et al. (1997a). 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 region of an enhanced number of particles (clump) 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, that often appear in the literature. While the former is normally associated with the fast kinetic response of resonant particles (typically of order ) studied in this work, the latter usually corresponds to a slow evolution (typically ) of Alfvénic modes, in the presence of a non-monotonic safety factor profile, that relies on the time variation of . 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 Sharapov et al. (2002). Chirping is faster and harder to suppress using external control.

Chirping modes can have frequency shifts greater than the linear growth rate and are observed to be precursors of even worst scenarios, known as avalanches. A spectrogram showing repetitive chirping cycles followed by avalanches for toroidal Alfvén eigenmodes (TAEs) in 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 Lilley et al. (2010). Fig. 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 bounce frequency to the power of Berk et al. (1997a). It has been successfully used for applications that include the inference of mode amplitude on MAST Pinches et al. (2004) and the estimation of kinetic parameters such as drive and damping in JT-60U Lesur et al. (2010) and in NSTX Fredrickson et al. (2006).

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, 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 we show how it is possible to include micro-turbulent stochasticity, which is shown to compete, and even greatly exceed, collisional scattering in many tokamak experiments and therefore needs to be added to the stochasticity introduced by pitch-angle scattering.

Chirping and quasilinear (QL) regimes correspond to two opposite 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 Drummond and Pines (1962); Vedenov et al. (1961) relies on averages, over a statistical ensemble, that 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 to be invoked. The fast-varying response associated to the ballistic term is disregarded, which imply 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 also is 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 Gorelenkov et al. (2016). It uses the usual structure of the QL system written in action-angle variables Kaufman (1972) with a broadened resonance width that scales with bounce frequency, growth rate and collisional frequency Berk et al. (1995); Ghantous et al. (2014).

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 is especially important for the development of the reduced models since their methodologies critically depend on that.

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 experimental analysis of 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.

## Ii Theoretical framework

### ii.1 Formulation of fast ion interaction with Alfvén waves

In axisymmetric tokamaks, , and (corresponding respectively to energy, canonical toroidal angular momentum and magnetic moment) 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 and are the mass and charge of EPs, is the EP speed, is the tokamak major radius, is the poloidal flux divided by , is the magnitude of the magnetic field, and are the poloidal and toroidal angles.

If the wave-particle interaction Hamiltonian is assumed, like in the NOVA code Gorelenkov et al. (1999a), to have dependences on time and poloidal and toroidal angle as , where and are the wave numbers associated to them, it follows that upon particle interaction with low-frequency modes, the invariances of and are broken but a new invariant arises: while is kept nearly constant. The existence of two invariants implies that the resonant particle dynamics is essentially one-dimensional Berk et al. (1997b). Therefore, a new variable (at constant ) can be used to describe this relevant one-dimensional path for EP dynamics (for steepest distribution function modification). The projection of the gradient operator onto this path is then given by

(1) |

The resonances for the several harmonics, which are multiple surfaces in space, are defined by

where is an integer and and are local toroidal and poloidal transit frequencies, respectively. The phase-space integration is represented as

where is the light speed and accounts for counter- and co-passing particles.

### ii.2 Collisional operator for fast ions

Collisional processes are an important element in the determination of the nonlinear character of wave oscillations Berk et al. (1996, 1997b); Lilley et al. (2009). 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 enhance 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 Goldston et al. (1981),

where is the distribution function; and are the pitch angle scattering rate and the inverse slowing down time, which are given by

(2) |

(3) |

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

(4) |

and

For a given ion species , , with being the proton mass, is the Coulomb logarithm, and are the background and energetic ion atomic number, is the density, is the temperature, is the mass and the thermal speed for a given species a is defined as . The three-dimensional velocity derivatives are projected onto the aforementioned path of the steepest gradient of the distribution function so that

Since and , 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 and drag coefficients are

(5) |

and

(6) |

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

### ii.3 Inclusion of micro-turbulent stochasticity

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 is often dominant in the determination of the global heat outflow. Micro-turbulence is introduced using a procedure that follows the one introduced by Lang and Fu Lang and Fu (2011), where it is 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 . Since with , the spatial micro-turbulence diffusion operator along the radial coordinate can be written as

(7) |

where is the EP diffusivity. Therefore, the ratio between the effective stochasticities coming from micro-turbulent (Eq. 7) and collisional processes (Eq. 5) is

(8) |

A subtlety in applying the model relies on the determination of . It is expected to be lower than thermal ion diffusivity 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 ions turbulent transport is negligible compared to Alfvénic-induced transport (see, for example, studies on ASDEX-U Geiger et al. (2015) and DIII-D Pace et al. (2013)), turbulence can be an important transport mechanism, as compared to collisions, at the onset of the evolution of modes, when their amplitude are still small. Over the past decade, the modeling of has been studied by several groups Hauff et al. (2009); Pueschel et al. (2012); Estrada-Mila et al. (2006); Albergante et al. (2010, 2009); Zhang et al. (2008). We have chosen to determine from the scalings that follow from gyrokinetic simulations of GTC Zhang et al. (2008), in which is proportional to the thermal ion diffusivity, , and a function of , where 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 to the inferred value of . Alternative ways of obtaining could be achieved by using Refs. Pueschel et al. (2012); Hauff et al. (2009). In our analysis we infer the value of from the outputs of the global transport code TRANSP Hawryluk (1980); Goldston et al. (1981) 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 to it because TRANSP cannot resolve well particle sources, especially close to wall. On the other hand, the thermal conductivity is reasonably well known and therefore is frequently used Heidbrink et al. (2009) as an indication of the actual value of (the exact relation for a Maxwellian distribution would be ).

### ii.4 The early nonlinear evolution of a mode amplitude

The onset of a mode amplitude evolution can be studied using perturbation theory within the kinetic framework. Refs. Berk et al. (1996); Breizman et al. (1997); Berk et al. (1997b); Lilley et al. (2009) showed that, for , truncation of mode amplitude at third order is justified due to fact that the unstable system is close to marginal stability. Taking (the overall stochasticity felt by EPs, which includes ) and independent of time but dependent on phase space localization, the equation for the early-time perturbed mode (a mode that exists without accounting for the kinetic component) amplitude evolution can be written as a time-delayed, integro-differential cubic equation

(9) |

where and

or

(10) |

accounts for the wave-particle energy exchange, where is the electric field eigenstructure and is the velocity of a resonant particle. are integers and are angles conjugated to the invariants of motion (actions of the Hamiltonian). In equation (9), the circumflex denotes normalization with respect to (growth rate minus damping rate) and is the time normalized with the same quantity. is the normalized complex mode amplitude of an eigenmode oscillating with frequency .

The solutions of eq. (9) can exhibit several bifurcations and therefore several phases, as shown for a bump-on-tail configuration in Ref. Berk et al. (1996). Interestingly, eq. (9) allows for the excitation of sub-critical instabilities and for nonlinear frequency splitting Fasoli et al. (1998). If the nonlinearity in equation (9) is weak, the system will most likely saturate close to the linear stability state, where the trapping frequency satisfies . 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 as the precursor to the formation of holes and clumps structures Berk et al. (1999). Furthermore, chirping events are significantly enhanced by the coherence introduced by dynamical friction (i.e, particle drag) Lilley et al. (2009, 2010) and are inhibited by stochasticity from diffusive processes, such as resonant particle heating Maslovsky et al. (2003), collisional pitch-angle scattering and from background turbulence, all of which contribute to causing particles to detune from a resonance. Stochastic events lead to loss of phase information that contribute to destroy coherent structures.

Eq. (9) was originally derived for a bump-on-tail system with Krook collisions Berk et al. (1996) and later generalized to complex tokamak geometries and also to include collisional scattering Berk et al. (1997b). Lilley, Breizman and Sharapov Lilley et al. (2009) 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.

### ii.5 A criterion for chirping onset

It has been shown Duarte et al. (2016) 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 would be that the real component of the right-hand side of Eq. (9) be negative at late times when the response is stationary, i.e. when the nonlinear term is allowed to balance the linear growth. The delta function that restricts the integration to the resonance condition can be exploited and the following criterion for the existence of fixed-frequency oscillations is obtained Duarte et al. (2016):

(11) |

where

(12) |

and is a normaliation for , which consists in the same sums and integrations that appear in the numerator of (11) but in the absence of . In eqs. (11) and (12) , the quantities , , , and are understood to be evaluated at . Criterion (11) was incorporated into NOVA-K making use of a polynomial interpolation for . 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 is plotted in figure (2) for different values of .

Drag introduces an oscillatory behavior of the integrand of . This has an effect of flipping the sign of the integral kernel of the cubic equation. This causes to vary more abruptly and therefore prevents a steady state solution from being achieved. Instead nonlinear chirping solution is more likely to be achieved. It is interesting to note that formally, drag enters the kinetic equation in a mathematically similar way as a chirping frequency does, with replaced by .

, which is a function of phase space, is plotted in Fig. (3) as a function of . Note that the positive and negative domains of are roughly an order of magnitude different in the interval , while for , .

For our investigation, we focus on the criterion for the non-existence of a fixed-frequency, steady solution. To determine the physical parameters needed in the calculation described below, the global transport code TRANSP Hawryluk (1980); Goldston et al. (1981) 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 due to 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 . 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.

## Iii Numerical study of the chirping criterion

### iii.1 Kinetic-MHD perturbative computations

NOVA Gorelenkov et al. (1999a) is a nonvariational, hybrid kinetic-MHD stability 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 Cheng (1992); Gorelenkov et al. (1999b) takes into account finite Larmor radius and orbit width efects to study the destabilization of MHD modes from the EPs free energy. We use NOVA-K to calculate the resonance surfaces, in 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 . 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 , with where is the mode eigenfrequency and is its growth rate, it is obtained Gorelenkov et al. (1999a); Cheng (1992)

where is the inertial (kinetic) energy, is the potential energy associated to the nonadiabatic component of the distribution function and is the power released by the resonant particles. The growth rate can be expressed as

(13) |

with the diamagnetic frequency being defined by

represents mode structure matrix elements, which are associated to the projection of the resonant particle current onto the wave electric field. For each isolated resonance, the reduced Hamiltionian can be written as Berk et al. (1997b)

where and the bounce frequency of the most deeply trapped particles can be shown to be

where the subscript denotes the resonance location and the derivative is defined by eq. (1). The matrices are related to (given by equation (10)) via

### iii.2 Mode structure identification

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 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 Zeeland et al. (2006). 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 Fredrickson et al. (2006); Podestà et al. (2012); Heidbrink (1995). This gives credence to the use of a linear code for the eigenstructure identification at early times.

### iii.3 Averaging implementations in NOVA-K

In order to calculate the above expressions for the effective collisional coefficients (5) and (6), 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.

#### iii.3.1 Bounce averaging

The period of a particle poloidal bounce or transit motion is given by , where the drift velocity is due to curvature and terms. Then, the bounce average of the effective collisional frequencies is given by

For each set , the trajectories are a priori known, hence there is no additional need to explicitly follow the particles motion within the code.

#### iii.3.2 Phase space averaging

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 to their relative contribution to the overall growth rate. Specifically, we evaluate

where is the contribution to the growth rate, , from a given phase space location that satifies . 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 Gorelenkov et al. (1999b).

### iii.4 Study of the chirping prediction for eigenmodes

In order to show the importance of micro-turbulence as a chirping suppression mechanism, an TAE mode 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 NOVA code is shown in Fig. (5). Shown in Fig. (6) are scans over the inferred experimental values of and , denoted by and , 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 and (as shown in Fig. (6)). The vertical axis represents the constants that are multiplying while the horizontal one represents constants that are multiplying . The point corresponds to the inferred experimental condition and the criterion boundary is represented by the countour. 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 sufficiently far from this transition boundary, steady solutions should arise in experiment.

Each contour plot corresponds to a given value of (eq. (11)), as labeled in the two plots in Fig. (6). In part (a), where is evaluated only in terms of collisional pitch-angle scattering and drag, is a negative number. Upon the addition of micro-turbulent scattering, the re-evaluation of shows that the point indicating the experimental conditions now operates where , 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 based on the scaling of Zhang et al. (2008) led to , which is consistent with direct measurements using beam blips in similar experiments carried out in TFTR Heidbrink and Sadler (1994). Numerically calculated and were and , respectively. was found to be nearly times larger than for this particular mode. The degree pitch-angle scattering and the inverse slowing down time were found to be and .

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 only depends on the ratio of the overall stochasticity and collisional drag. From equations (2), (3), (5) and (6), we see that the effective drag coefficient is proportional to while the effective collisional scattering is insensitive to . It is then possible to mimic the inclusion of turbulence into by artificially decreasing drag, which can be achieved by multiplying by a corresponding numerical factor that leads to the correct value of the averaged values of .

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 (), and strong toroidal rotation (up to 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) Gorelenkov et al. (2007) 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 ( 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

The time-delayed cubic equation (9) is derived assuming that the nonlinear bounce frequency is much less than the growth rate in absence of dissipation. 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 , is still small compared to . From DIII-D observations, we have noted that the chirping behavior typically starts when drops to values lower than Duarte et al. (2016) but chirping persists even when 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 .

Chirping does not restart when drops for the second time, at . 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 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.

Fig (10) (a) shows the experimental condition for the mode in DIII-D shot 152828 shown in Fig. (9) before chirping starts (at , when ) and Fig. (10) (b) shows it during chirping (at , when ), 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 , used the mode structure calculated at (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 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 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 bucket structures that sustain chirping.

## Iv Summary and discussions

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 for the general problem, with a similar perturbative procedure used in Berk et al. (1997b). 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 be 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. (9)).

An interesting effect unveiled by the methodology used in this paper is that the disparity of magnitudes of the positive and negative regions of (equation (12)), typically leads to the prediction that a given mode should chirp, even for , unless a strong diffusive mechanism implies . 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 is larger for . 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 then 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.

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. Zhang et al. (2008), based on electrostatic turbulence. However, the value chosen is rough, and the actual scaling for a given experiment might differ considerably from our choice. Using the estimates for EP diffusivity, we brought in the model the effect of background turbulence, that adds up to collisional scattering and contribute to suppress chirping.

A number of factors that may influence the chirping formation are not captured by the current theory. For example, static 3D fields Bortolon et al. (2013) 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 Maslovsky et al. (2003); Heidbrink et al. (2006), 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.

## Appendix A Chirping likelihood in terms of the beam injection speed

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 and , 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, while fixed-frequency, quasi-steady modes are in general observed when the NBI is sub-Alfvénic Kusama et al. (1999); Kramer et al. (2000); Shinohara et al. (2002, 2004); Ishikawa et al. (2005). The ALEs have been intensively studied over the past decade Briguglio et al. (2007); Bierwage et al. (2013, 2014, 2016); Bierwage and Shinohara (2016); Bierwage et al. (2017). 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 a TAE for which the resonance condition is . For scenarios in which micro-turbulence can be ignored in comparison with collisional scattering, the relevant ratio for the criterion, eq. (11), would be simply , with the effective scattering given by eq. (5),

where is the EP cyclotron frequency and . The effective drag is given by (6),

with given by (4). Therefore

(14) |

where

and is the angle between and . The ratio between the chirping-criterion-relevant parameter prescribed by eq. (14), for the cases in which the main resonance of a TAE is at and is plotted in Fig. (13). The ratio is for . If , it means that supra-Alfvénic injection implies more likelihood for chirping. It is a direct implication of the chirping criterion (eq. (11)) and the form of (Fig. (12)). Considering deuterium as the only background ion species, , , and using typical parameters for JT-60U, such as , and one obtains . The factor contributes to reduce the variable 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. 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 Zhang et al. (2008) does not change the conclusion that higher resonant velocities imply more likelihood for chirping, provided that is sufficiently small. In fact, turbulence even contributes to strenghten this interpretation. This is because Zhang et al. (2008) found that the micro-turbulence diffusivity goes as for passing particles and 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 , which makes chirping for supra-Alfvénic injection even more probable than in the purely collisional case.

###### Acknowledgements.

We acknowledge fruitful discussions with G.-Y. Fu, E. D. Fredrickson and R. M. O. Galvão. This work was supported by the São Paulo Research Foundation (FAPESP, Brazil) under grants 2012/22830-2 and 2014/03289-4, and by US Department of Energy (DOE) under contracts 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”.## References

- Heidbrink (2008) W. W. Heidbrink, Phys. Plasmas 15, 055501 (2008), http://dx.doi.org/10.1063/1.2838239.
- Gorelenkov et al. (2014) N. Gorelenkov, S. Pinches, and K. Toi, Nucl. Fusion 54, 125001 (2014).
- Lauber (2013) P. Lauber, Physics Reports 533, 33 (2013).
- Chen and Zonca (2016) L. Chen and F. Zonca, Rev. Mod. Phys. 88, 015008 (2016).
- Berk et al. (1996) H. L. Berk, B. N. Breizman, and M. Pekker, Phys. Rev. Lett. 76, 1256 (1996).
- Berk et al. (1997a) H. Berk, B. Breizman, and N. Petviashvili, Phys. Lett. A 234, 213 (1997a).
- Bernstein et al. (1957) I. B. Bernstein, J. M. Greene, and M. D. Kruskal, Phys. Rev. 108, 546 (1957).
- Sharapov et al. (2002) S. E. Sharapov, B. Alper, H. L. Berk, D. N. Borba, B. N. Breizman, C. D. Challis, A. Fasoli, N. C. Hawkes, T. C. Hender, J. Mailloux, S. D. Pinches, D. Testa, and E.-J. work programme, Physics of Plasmas 9, 2027 (2002).
- Lilley et al. (2010) M. K. Lilley, B. N. Breizman, and S. E. Sharapov, Phys. Plasmas 17, 092305 (2010), http://dx.doi.org/10.1063/1.3486535.
- Pinches et al. (2004) S. D. Pinches, H. L. Berk, M. P. Gryaznevich, S. E. Sharapov, and J.-E. Contributors, Plasma Phys. Control. Fusion 46, S47 (2004).
- Lesur et al. (2010) M. Lesur, Y. Idomura, K. Shinohara, X. Garbet, and the JT-60 Team, Phys. Plasmas 17, 122311 (2010), http://dx.doi.org/10.1063/1.3500224.
- Fredrickson et al. (2006) E. D. Fredrickson, R. E. Bell, D. S. Darrow, G. Y. Fu, N. N. Gorelenkov, B. P. LeBlanc, S. S. Medley, J. E. Menard, H. Park, A. L. Roquemore, W. W. Heidbrink, S. A. Sabbagh, D. Stutman, K. Tritz, N. A. Crocker, S. Kubota, W. Peebles, K. C. Lee, and F. M. Levinton, Phys. Plasmas 13, 056109 (2006), http://dx.doi.org/10.1063/1.2178788.
- Drummond and Pines (1962) W. Drummond and D. Pines, Nucl. Fusion Suppl. 2, Pt. 3 (1962).
- Vedenov et al. (1961) A. A. Vedenov, E. P. Velikhov, and R. Z. Sagdeev, Sov. Phys. Uspekhi 4, 332 (1961).
- Gorelenkov et al. (2016) N. Gorelenkov, V. Duarte, and H. Berk, in APS Meeting Abstracts (2016).
- Kaufman (1972) A. N. Kaufman, Phys. Fluids 15, 1063 (1972).
- Berk et al. (1995) H. Berk, B. Breizman, J. Fitzpatrick, and H. Wong, Nucl. Fusion 35, 1661 (1995).
- Ghantous et al. (2014) K. Ghantous, H. L. Berk, and N. N. Gorelenkov, Phys. Plasmas 21, 032119 (2014), http://dx.doi.org/10.1063/1.4869242.
- Gorelenkov et al. (1999a) N. N. Gorelenkov, C. Z. Cheng, and G. Y. Fu, Phys. Plasmas 6, 2802 (1999a).
- Berk et al. (1997b) H. L. Berk, B. N. Breizman, and M. Pekker, Plasma Phys. Rep. 23, 778 (1997b).
- Lilley et al. (2009) M. K. Lilley, B. N. Breizman, and S. E. Sharapov, Phys. Rev. Lett. 102 (2009), 10.1103/physrevlett.102.195003.
- Goldston et al. (1981) R. J. Goldston, D. C. McCune, H. H. Towner, L. S. Davis, R. J. Hawryluk, and G. L. Schmidt, J. Comput. Phys. 43, 61 (1981).
- Lang and Fu (2011) J. Lang and G.-Y. Fu, Phys. Plasmas 18, 055902 (2011).
- Geiger et al. (2015) B. Geiger, M. Weiland, A. Mlynek, M. Reich, A. Bock, M. Dunne, R. Dux, E. Fable, R. Fischer, M. Garcia-Munoz, J. Hobirk, C. Hopf, S. Nielsen, T. Odstrcil, C. Rapson, D. Rittich, F. Ryter, M. Salewski, P. A. Schneider, G. Tardini, and M. Willensdorfer, Plasma Phys. Control. Fusion 57, 014018 (2015).
- Pace et al. (2013) D. C. Pace, M. E. Austin, E. M. Bass, R. V. Budny, W. W. Heidbrink, J. C. Hillesheim, C. T. Holcomb, M. Gorelenkova, B. A. Grierson, D. C. McCune, G. R. McKee, C. M. Muscatello, J. M. Park, C. C. Petty, T. L. Rhodes, G. M. Staebler, T. Suzuki, M. A. Van Zeeland, R. E. Waltz, G. Wang, A. E. White, Z. Yan, X. Yuan, and Y. B. Zhu, Phys. Plasmas 20, 056108 (2013).
- Hauff et al. (2009) T. Hauff, M. J. Pueschel, T. Dannert, and F. Jenko, Phys. Rev. Lett. 102, 075004 (2009).
- Pueschel et al. (2012) M. Pueschel, F. Jenko, M. Schneller, T. Hauff, S. Guenter, and G. Tardini, Nuclear Fusion 52, 103018 (2012).
- Estrada-Mila et al. (2006) C. Estrada-Mila, J. Candy, and R. E. Waltz, Physics of Plasmas 13, 112303 (2006).
- Albergante et al. (2010) M. Albergante, J. Graves, A. Fasoli, and X. Lapillonne, Nuclear Fusion 50, 084013 (2010).
- Albergante et al. (2009) M. Albergante, J. P. Graves, A. Fasoli, F. Jenko, and T. Dannert, Physics of Plasmas 16, 112301 (2009), http://dx.doi.org/10.1063/1.3257913 .
- Zhang et al. (2008) W. Zhang, Z. Lin, and L. Chen, Phys. Rev. Lett. 101 (2008), 10.1103/physrevlett.101.095001.
- Hawryluk (1980) R. J. Hawryluk, in Physics of Plasmas Close to Thermonuclear Conditions, Vol. 1, edited by B. Coppi, G. G. Leotta, D. Pfirsch, R. Pozzoli, and E. Sindoni (CEC, Brussels, 1980) pp. 19–46.
- Heidbrink et al. (2009) W. W. Heidbrink, J. M. Park, M. Murakami, C. C. Petty, C. Holcomb, and M. A. V. Zeeland, Phys. Rev. Lett. 103 (2009), 10.1103/physrevlett.103.175001.
- Breizman et al. (1997) B. N. Breizman, H. L. Berk, M. S. Pekker, F. Porcelli, G. V. Stupakov, and K. L. Wong, Phys. Plasmas 4, 1559 (1997).
- Fasoli et al. (1998) A. Fasoli, B. N. Breizman, D. Borba, R. F. Heeter, M. S. Pekker, and S. E. Sharapov, Phys. Rev. Lett. 81, 5564 (1998).
- Berk et al. (1999) H. L. Berk, B. N. Breizman, J. Candy, M. Pekker, and N. V. Petviashvili, Phys. Plasmas 6, 3102 (1999).
- Maslovsky et al. (2003) D. Maslovsky, B. Levitt, and M. E. Mauel, Phys. Rev. Lett. 90, 185001 (2003).
- Duarte et al. (2016) V. Duarte, H. Berk, N. Gorelenkov, W. Heidbrink, G. Kramer, R. Nazikian, D. Pace, M. Podesta, B. Tobias, and M. V. Zeeland, “Prediction of nonlinear evolution character of energetic-particle-driven instabilities,” (2016), arXiv:1610.05059 .
- Cheng (1992) C. Cheng, Phys. Rep. 211, 1 (1992).
- Gorelenkov et al. (1999b) N. N. Gorelenkov, Y. Chen, R. B. White, and H. L. Berk, Phys. Plasmas 6, 629 (1999b).
- Zeeland et al. (2006) M. A. V. Zeeland, G. J. Kramer, M. E. Austin, R. L. Boivin, W. W. Heidbrink, M. A. Makowski, G. R. McKee, R. Nazikian, W. M. Solomon, and G. Wang, Phys. Rev. Lett. 97 (2006), 10.1103/physrevlett.97.135001.
- Podestà et al. (2012) M. Podestà, R. Bell, A. Bortolon, N. Crocker, D. Darrow, A. Diallo, E. Fredrickson, G.-Y. Fu, N. Gorelenkov, W. Heidbrink, G. Kramer, S. Kubota, B. LeBlanc, S. Medley, and H. Yuh, Nucl. Fusion 52, 094001 (2012).
- Heidbrink (1995) W. W. Heidbrink, Plasma Phys. Control. Fusion 37, 937 (1995).
- Heidbrink and Sadler (1994) W. W. Heidbrink and G. J. Sadler, Nuclear Fusion 34, 535 (1994).
- Fu et al. (1998) G. Y. Fu, R. Nazikian, R. Budny, and Z. Chang, Physics of Plasmas 5, 4284 (1998), http://aip.scitation.org/doi/pdf/10.1063/1.873165 .
- Gorelenkov et al. (2007) N. Gorelenkov, H. Berk, E. Fredrickson, S. Sharapov, and J. E. Contributors, Physics Letters A 370, 70 (2007).
- Bortolon et al. (2013) A. Bortolon, W. W. Heidbrink, G. J. Kramer, J.-K. Park, E. D. Fredrickson, J. D. Lore, and M. Podestà, Phys. Rev. Lett. 110, 265008 (2013).
- Heidbrink et al. (2006) W. W. Heidbrink, E. Ruskov, E. D. Fredrickson, N. Gorelenkov, S. S. Medley, H. L. Berk, and R. W. Harvey, Plasma Physics and Controlled Fusion 48, 1347 (2006).
- Kusama et al. (1999) Y. Kusama, G. Kramer, H. Kimura, M. Saigusa, T. Ozeki, K. Tobita, T. Oikawa, K. Shinohara, T. Kondoh, M. Moriyama, F. Tchernychev, M. Nemoto, A. Morioka, M. Iwase, N. Isei, T. Fujita, S. Takeji, M. Kuriyama, R. Nazikian, G. Fu, K. Hill, and C. Cheng, Nucl. Fusion 39, 1837 (1999).
- Kramer et al. (2000) G. Kramer, M. Iwase, Y. Kusama, A. Morioka, M. Nemoto, T. Nishitani, K. Shinohara, S. Takeji, K. Tobita, T. Ozeki, C. Cheng, G.-Y. Fu, and R. Nazikian, Nucl. Fusion 40, 1383 (2000).
- Shinohara et al. (2002) K. Shinohara, M. Takechi, M. Ishikawa, Y. Kusama, A. Morioka, N. Oyama, K. Tobita, T. Ozeki, the JT-60 Team, N. Gorelenkov, C. Cheng, G. Kramer, and R. Nazikian, Nuclear Fusion 42, 942 (2002).
- Shinohara et al. (2004) K. Shinohara, M. Takechi, M. Ishikawa, Y. Kusama, K. Tsuzuki, K. Urata, H. Kawashima, K. Tobita, A. Fukuyama, C. Z. Cheng, D. S. Darrow, G. J. Kramer, N. N. Gorelenkov, R. Nazikian, Y. Todo, Y. Miura, and T. Ozeki, Plasma Physics and Controlled Fusion 46, S31 (2004).
- Ishikawa et al. (2005) M. Ishikawa, M. Takechi, K. Shinohara, Y. Kusama, C. Cheng, G. Matsunaga, Y. Todo, N. Gorelenkov, G. Kramer, R. Nazikian, A. Fukuyama, V. Krasilnikov, Y. Kashuck, T. Nishitani, A. Morioka, M. Sasao, M. Isobe, and the JT-60 team, Nuclear Fusion 45, 1474 (2005).
- Briguglio et al. (2007) S. Briguglio, G. Fogaccia, G. Vlad, F. Zonca, K. Shinohara, M. Ishikawa, and M. Takechi, Physics of Plasmas 14, 055904 (2007), http://dx.doi.org/10.1063/1.2710208 .
- Bierwage et al. (2013) A. Bierwage, K. Shinohara, N. Aiba, and Y. Todo, Nuclear Fusion 53, 073007 (2013).
- Bierwage et al. (2014) A. Bierwage, Y. Todo, N. Aiba, and K. Shinohara, Nuclear Fusion 54, 104001 (2014).
- Bierwage et al. (2016) A. Bierwage, Y. Todo, N. Aiba, and K. Shinohara, Nuclear Fusion 56, 106009 (2016).
- Bierwage and Shinohara (2016) A. Bierwage and K. Shinohara, Physics of Plasmas 23, 042512 (2016), http://dx.doi.org/10.1063/1.4947034 .
- Bierwage et al. (2017) A. Bierwage, K. Shinohara, Y. Todo, N. Aiba, M. Ishikawa, G. Matsunaga, M. Takechi, and M. Yagi, Nuclear Fusion 57, 016036 (2017).