High harmonic ion cyclotron heating in DIII-D: Beam ion absorption and sawtooth stabilization

Combined neutral beam injection and fast wave heating at the fourth cyclotron harmonic produce an energetic deuterium beam ion tail in the DIII-D tokamak. When the concentration of thermal hydrogen exceeds ∼5%, the beam ion absorption is suppressed in favour of second harmonic hydrogen absorption. As theoretically expected, the beam absorption increases with beam ion gyro-radius; also, central absorption at the fifth harmonic is weaker than central absorption at the fourth harmonic. For central heating at the fourth harmonic, an energetic, perpendicular, beam population forms inside the q = 1 surface. The beam ion tail transiently stabilizes the sawtooth instability but destabilizes toroidicity induced Alfvén eigenmodes (TAEs). Saturation of the central heating correlates with the onset of the TAEs. Continued expansion of the q = 1 radius eventually precipitates a sawtooth crash; complete magnetic reconnection is observed.


Introduction
Cyclotron damping of fast waves in the ion cyclotron range of frequencies is a standard heating scheme in tokamaks. Absorption at the fundamental cyclotron frequency by a low concentration 'minority' species has formed an energetic, anisotropic 'tail' distribution in many tokamaks [1]. Tails produced by second harmonic heating are also well documented [1] but less information is available about high harmonic heating. Second, third and fourth harmonic acceleration of hydrogen were studied on Tokapole II [2]. At JET, strong third harmonic acceleration of (initially) thermal deuterium ions was observed [3]. Beam ion acceleration during combined beam injection and third harmonic ICRF heating was measured on JIPP TII-U [4], JT-60 [5] and TEXTOR [6,7], while fourth harmonic heating of beam ions was briefly studied on JT-60 [5] and JET [8]. Hydrogen beam ions were accelerated by second, third and fourth harmonic heating at JT-60U, but second harmonic heating was most efficient [9]. Beam ion absorption at high ion cyclotron harmonics was also observed during multi-frequency fast wave current drive (FWCD) experiments on DIII-D [10]; the beam ion absorption was originally attributed to sixth and seventh harmonic damping at 83 MHz, but the results presented here indicate that fourth harmonic damping at 60 MHz is a more likely explanation.
This article discusses fourth harmonic cyclotron heating of deuterium beam ions in the DIII-D tokamak and the effect of the accelerated beam population on sawtooth and TAE stability. The documentation of fourth harmonic heating in a tokamak with good tail ion confinement is the most extensive yet reported. The study of the transiently stabilized sawtooth cycle is the first to rely on q profiles from an absolutely calibrated motional Stark effect (MSE) diagnostic.
The following interpretive scenario is consistent with the data. With the cyclotron resonance near the magnetic axis, fourth harmonic absorption on the deuterium beam ions dominates over second harmonic absorption on the ambient (º2%) thermal protons, primarily because of the relatively large gyro-radius of the beam ions (Section 2). An energetic, perpendicular, beam population forms inside the q = 1 surface (Section 2). The precessing beam ions transiently stabilize the internal kink (Section 3.1) but, when the beam population is sufficiently large, TAEs are excited (Section 3.2). The TAEs limit the central electron temperature and, presumably, the central beam density. Meanwhile, current accumulates on-axis and the q = 1 radius expands, until kink stability is lost (Section 3.3). To eliminate sawtooth crashes altogether, current profile control is required (Section 4).

Experimental observations
The data in this article are from several different days of the 1998 experimental campaign. During 1998, the DIII-D tokamak [16] (major radius R 0 1.7 m, minor radius a 0.6 m, graphite walls, deuterium plasmas) was equipped with three four-element antenna arrays that were connected to three 2 MW generators [17]. For most of these experiments, one transmitter operated at ∼60 MHz while the other two operated at ∼83 MHz; typical coupled power levels into L mode plasmas were 1 MW at each frequency. For most of the discharges, the outer gap (the distance from the last closed flux surface to the limiter at the outboard midplane) was 3-4 cm, although outer gap spacings as large as 9 cm were occasionally used. For the 60 MHz system, countercurrent drive phasing (90 • toroidal phasing between straps) was normally employed. The peak in the vacuum spectrum was at n = k c/ω 5. At the usual toroidal field of 1.9 T, the 60 MHz waves resonated with the deuterium cyclotron motion at several radial locations ( Fig. 1), with the central resonance corresponding to the fourth harmonic. Most plasmas were upper single null divertor discharges with an elongation κ 1.75. The ∇B drift was downwards (away from the X point) in an attempt to avoid H mode transitions. The beams injected 70-86 keV deuterium neutrals in the direction of the plasma current at tangency radii of 0.76 m (for the so-called 'Right' sources) and 1.15 m (for the 'Left' sources).
The temporal evolution of a typical discharge is shown in Fig. 2. Neutral beam injection commences 0.25 s after the current reaches its flat top value of 1.2 MA. To accommodate the diagnostics, five modulated sources inject into the plasma but the total average injected power of 1.6 MW is less than that injected by a single steady source. To minimize perturbations in the beam population, the combined beam waveform never turns off for more than 10 ms, which is much less than the beam ion slowing down time of τ s 75 ms. The 60 MHz system couples ∼1 MW of power between 1.8 and 3.0 s. With the application of RF power, the stored energy, d-d neutron emission [18] and flux of energetic perpendicular neutrals [19] increase on a τ s timescale. ECE measurements with a heterodyne radiometer [20] indicate that the central electron temperature increases during RF heating and that the sawtooth cycle is modified.
Thomson scattering [21] and interferometer [22] data show that the line averaged electron density is typicallyn e = 3 × 10 13 cm −3 in the L mode plasmas. Measurements with the charge exchange recombination (CER) diagnostic [23] during RF heating indicate that the central ion temperature is ∼2.6 keV, the central toroidal rotation frequency is ∼7 kHz and the central carbon density n C (0) 1.3 × 10 12 cm −3 , which implies that the effective charge is Z eff 2.0. A pair of magnetic loops (bandwidth > 1 MHz) that are mounted ∼0.8 m above the outer midplane and separated by 15 • toroidally monitors Alfvén activity with a frequency response up to 1 MHz.
Equilibria are reconstructed using the EFIT code [24]. For this study, the 'standard' reconstructions fit data from the absolutely calibrated, 36 channel MSE diagnostic [25] and the arrays of external magnetic probes.
In discharge 96033 (Fig. 2), the average beam injection energy E is 74 keV and the effect of the RF heating on the plasma parameters is relatively modest. For a similar discharge (No. 96043) with E = 80 keV, the increases in perpendicular stored energy, neutron rate and sawtooth period are much greater, as will be shown.  . Time evolution of the measured neutron emission (solid curve) and the classically expected rate as predicted by TRANSP [26] (dotted curve) and by a zero dimensional code [18] (dashed curve) in discharge 96043.
As expected for cyclotron damping, the perpendicular stored energy increases faster than the parallel stored energy during the RF pulse (Fig. 3). (The diamagnetic loop signal was imperfectly compensated during the 1998 campaign, but relative changes during the current flat top correlate with other measurements of the fast ion population.) The plasma anisotropy ∆β ⊥ = β dia − β eq is approximately 10% of the plasma beta during RF heating.
Neutron measurements show that deuterium ions are accelerated. Predictions of the d(d, n) rate that assume classical beam ion deposition and thermalization are consistent with the measured rate prior to the RF pulse (Fig. 4). During the RF heating, however, the neutron rate is strongly modulated . Electron temperature profile after mapping to flux co-ordinates for the kinetic EFIT reconstruction at (a) 1.68 s and (b) 2.00 s in discharge 96043. ECE data from outside the magnetic axis (open diamonds), inside the magnetic axis (solid diamonds) and data from Thomson scattering (crosses) are shown. During the RF pulse, the central pressure is underestimated, resulting in misalignment of the ECE data. (c) The measured thermal pressure (solid curve) and classical beam pressure (dotted curve) at 2.00 s, and an adjusted beam pressure (dashed curve) that yields properly aligned ECE data. The error bar represents the random error associated with uncertainties in ne (the dominant error), Te, Ti and nC . The abscissa is the normalized square root of the toroidal flux ρ. by the sawtooth cycle, increasing to values 40-50% larger than the classical rate prior to the sawtooth crash. The difference between the measured and predicted signals is due to the neglect of ion cyclotron acceleration in the calculations. Under these conditions, TRANSP (Ref. [26] and references therein) calculations indicate that beam-plasma reactions classically account for ∼90% of the reaction rate (beambeam reactions contribute ∼10% and the thermonuclear rate is º1%), so the enhanced neutron rate is caused by acceleration of the beam ions.
Detailed analysis of the equilibrium [10] shows that the pressure increases near the magnetic axis. In the 'standard' EFIT reconstruction, p and f f are adjusted to minimize the χ 2 value of the magnetics and MSE data. (Here, p is the pressure gradient, f is the poloidal current enclosed between a flux surface and the central symmetric axis and χ 2 is the usual statistical deviation between measured and modelled quantities.) In a 'kinetic' EFIT reconstruction, measurements of the pressure profile are also included in the minimization. The thermal pressure profile is obtained from measurements of T e , T i , n e and n C , but the fast ion pressure profile is not measured. The usual procedure is to assume classical beam deposition and thermalization to obtain a calculated beam pressure profile; however, this procedure neglects RF acceleration of the beam ions. The classical beam distribution yields an acceptable equilibrium reconstruction prior to the RF pulse but an incorrect equilibrium reconstruction during the RF (Fig. 5). To achieve acceptable agreement with the MSE data and the ECE measurement of the magnetic axis, the central beam pressure must be doubled or tripled relative to the classical value. Integration of the adjusted pressure profile indicates that the equilibrium stored energy (which weights the parallel and perpendicular pressures equally) is ²5% larger than classically expected during the RF pulse.
Active charge exchange measurements confirm that energetic, perpendicular, centrally located ions are accelerated. The geometry of the four charge exchange lines of sight are illustrated in Fig. 1. All four analysers measure trapped ions. (The three vertical analysers measure ions with no toroidal velocity, while the horizontal analyser has a tangency radius of ∼0.25 m for these experiments.) No energetic tail ions are detected by the outermost (V3) analyser at R = 2.1 m. Spatially resolved (active) measurements obtained by modulating a heating beam are available for the other three analysers (Fig. 1). All three of these analysers detect 30-60 keV neutrals during beam injection, but the flux from R 1.35 m (H) and R 1.5 m (V1) does not change measurably during the RF pulse. The only signal that changes appreciably is from the analyser that measures neutrals born predominately inside the q = 1 surface (the analyser V2 at R 1.94 m), confirming the central nature of the RF heating. An example of the unprocessed signal from the analyser V2 is shown in Fig. 2. The spikes in the neutral flux occur when one of the diagnostic heating beams is pulsed; this is the active charge exchange signal. Figure 6 shows the V2 energy spectrum obtained by averaging data from six similar L mode plasmas with average injection energies of 74-76 keV, average injected beam powers of 1.4-2.5 MW and coupled 60 MHz RF powers RF ON RF OFF ENERGY (keV) ACTIVE FLUX (V) Figure 6. The V2 neutral particle spectra from six discharges similar to discharge 96033 during beam heating alone (open diamonds) and during combined beam and 60 MHz RF heating (crosses). The active flux is the difference between the flux when the diagnostic heating beam is on and when it is off. The analyser V2 measures neutrals born just outside the magnetic axis with zero toroidal velocity (Fig. 1). The error bars represent statistical variations in the binned data. of 0.9-1.1 MW. The RF heating affects higher energy beam ions more strongly than lower energy ions.
[Corrections for attenuation of the diagnostic beam and escaping neutrals are neglected in Fig. 6 but, because the density is 10-25% higher during the RF pulse, these corrections increase the difference between the 'RF on' and 'RF off' points. (The calculated change is º50%.)] Presumably, beam ions are accelerated to energies well above the injection energy, but the analyser used in the 1998 campaign was unable to measure neutrals above 75 keV.
The ion cyclotron waves heat the bulk plasma effectively in this regime. The energy confinement time with beams alone (prior to the RF) averages 84 ms for discharges like No. 96033 (Fig. 2). With the addition of the RF, the incremental confinement time, τ inc ≡ ∆W/∆P = ∆W eq /(P RF + I P ∆V L ), averages 104 ms. The total radiated power increases by ∼300 kW with 1 MW of 60 MHz heating, but most of this is from the divertor; the core radiated power barely increases (º70 kW).
As well as fourth harmonic absorption by the deuterium beam ions, the fast waves can be absorbed by thermal hydrogen at the second harmonic. To investigate this effect, hydrogen gas was puffed into the vessel in a series of discharges similar to the one shown in Fig. 2. Absorption by beam ions results in an enhancement of the neutron rate S exp relative to the classically expected rate S cl (Fig. 4). As a convenient measure of beam ion damping, we define the neutron enhancementŜ as the ratio of the peak   [29] of power absorbed by beam ions (squares), thermal hydrogen (triangles) and electrons (crosses). Each point represents the absorption after three ray passes in the measured experimental profiles. measured rate to the zero dimensional classical prediction [18] during the RF heating, normalized to the same ratio prior to the RF heating, When hydrogen gas is puffed into the vessel, the neutron enhancement steadily decreases ( Fig. 7(a)). Only a modest rise in electron density (º10%) is associated with the gas puff, so the reduction in neutron rate is attributed to competitive wave absorption by the thermal hydrogen ( Fig. 7(c)). The central electron temperature follows the same trend, decreasing from 2.6 to 1.9 keV in the scan; in contrast, the ion temperature changes little ( Fig. 7(b)). The incremental confinement time τ inc decreases ∼25% at the largest hydrogen concentration. The reason for this reduction is not certain. One possibility is that thermal hydrogen is more poorly confined than beam ions in these L mode plasmas. [Studies in several tokamaks [1] find that the beam ion diffusivity is typically º0.1 m 2 /s, but the deuterium thermal diffusivity (derived from power balance analysis) is roughly 3 m 2 /s at ρ = 0.5 in these plasmas.] The shortening of the sawtooth period may also play a role. The efficiency of tail formation depends sensitively on antenna phasing. Generally the toroidal phase was chosen to drive current opposite to the plasma current (called counter-phasing at DIII-D); the waves (and associated electron drift) propagate in the direction of the plasma current. On one discharge, however, the phasing was reversed to drive positive current. Comparison of two virtually identical discharges (same shape and profiles of q, n e , T e and T i prior to the wave heating) shows that the discharge with counter-phasing has stronger central electron and beam ion heating than the discharge with co-phasing (Fig. 8). The difference in stored energy is about 4%.
The data in Figs 1-7 are all from similar discharges formed on the same day. To investigate the generality of these observations, a database of all the 60 MHz data from the 1998 campaign was assembled. Some data acquired at 83 MHz were also included. Figure 9 shows that the largest increases in neutron emission occur when the fourth harmonic resonance layer is situated near the magnetic axis. Fifth harmonic heating results in smaller tails.
A closer examination of the data in Fig. 9 indicates that the frequency of maximum absorption is affected by the Doppler shift. Restricting the data to the 36 discharges with an 80 keV tangential (Left) source and over 0.7 MW of 60 MHz power, the largest values of neutron enhancement occur when ω/Ω D (0) 4.14. The observed shift of 0.14 is in good agreement with the expected Doppler shift for damping on co-passing beam ions [27], ln v /c 0.16, so the maximum absorption occurs when the Doppler shifted resonance is near the magnetic axis, as expected.
The gyro-radius dependence of tail formation was studied by injecting higher energy, more perpendicular (Right) sources than the 80 keV, tangential (Left) sources used in most of the discharges. An increase in the injected gyro-radius by 25% results in ∼50% larger neutron enhancement factors (Fig. 10). This positive dependence of tail heating on beam ion gyroradius is evident for the entire database of L mode plasmas from the 1998 campaign (correlation coefficient [28] To summarize, the principal experimental results are: (e) Second harmonic heating of thermal hydrogen is significant for hydrogen concentrations of O(10%).

Theory of wave absorption
A complete theoretical treatment of the expected absorption in these conditions requires coupled solutions of the wave equations and of the beam ion Fokker-Planck equation. In future campaigns, we plan to measure the beam distribution in the energy range 0.1-1.0 MeV. Comparison of these measurements with full solutions of the Fokker-Planck equation will be reported in a future publication. As a simple alternative, we use ray tracing to study the wave physics and a local solution of the Fokker-Planck equation to study the beam ion distribution function. This simple approach successfully accounts for all of the experimental trends documented in Section 2.1.
For these experiments, the wavelength is much smaller than both the minor radius (λ/a 0.16) and the density gradient scale length in the plasma core, so ray tracing techniques are valid. The CUR-RAY code [29] represents the launched wave spectrum by an array of rays with power weighted phasing n , then uses a warm plasma dispersion relation for fast waves to follow the rays in the experimental equilibrium. Finite-beta correction terms [30] are included in the dispersion relation. Power absorption is calculated for electron Landau and TTMP damping, thermal ion cyclotron damping (all harmonics) and beam ion cyclotron damping (all harmonics) [31]. The beam population is approximated as an isotropic Maxwellian distribution, using the classical distribution found by the ONETWO transport code [32] in the absence of RF acceleration (weak RF limit). (Calculations with the anisotropic slowing down distribution computed by TRANSP [26] yield similar results.) The calculation of high harmonic ion cyclotron damping includes exact evaluation of the modified Bessel functions and harmonic overlap effects. Simple analytical formulas developed by Porkolab [31] provide a check on the CURRAY results.
The results of the analytical and ray tracing calculations for a typical L mode plasma (similar to the one shown in Fig. 2) are compiled in Tables 1 and 2. Two cases are compared: 60 MHz heating (fourth harmonic for the beam ions) and 83 MHz heating (fifth harmonic). The calculated ray trajectories are qualitatively similar in the two cases; in both cases,   central absorption dominates on the first few passes of the rays. More of the launched power is absorbed at 60 MHz than at 83 MHz, however (58% on the first few passes versus 19%). Both the simple analytical calculation and the ray tracing calculation find that beam ion absorption is dominant at 60 MHz, consistent with the experimental observations of strong fourth harmonic beam ion acceleration. Both calculations also successfully predict the weaker beam ion acceleration at 83 MHz (Tables 1 and 2). The predictions have a simple mathematical explanation: even for k ⊥ ρ i ∼ 1, high order Bessel functions decrease rapidly with increasing harmonic number. The ray tracing calculations also account for the results of the hydrogen gas puffing experiment. As the hydrogen concentration increases, the hydrogen absorption increases and both the beamion and electron absorption decrease (Fig. 7(c)), which explains the observed reduction in neutron rate ( Fig. 7(a)). The observed reduction in electron temperature ( Fig. 7(b)) is probably caused by the reduction in direct absorption and by a reduction in collisional power transfer from the tail ions. (Because of their high velocities, accelerated beam ions collide primarily with electrons, but accelerated thermal protons transfer energy to both electrons and ions.) The basic reason fourth harmonic absorption by beam ions is larger than second harmonic absorption by thermal hydrogen is that the beam ion gyroradius is nine times larger than the thermal hydrogen gyro-radius.
For the experimental case shown in Fig. 8, the ray tracing code predicts that 30% of the RF power is absorbed by beam ions with counter-phasing while, with co-phasing, the predicted beam ion absorption is only 22% (Table 3). In this calculation, the classical slowing down distribution computed by TRANSP [26] is employed. Both the ray trajectories and the location of the Doppler shifted resonance layer are affected by the phasing. For counter-phasing, the Doppler shift moves the resonance layer ∼7 cm closer to the magnetic axis, but for co-phasing the shift is in the opposite direction (Fig. 9). The prediction is consistent with the experimental trends, but the observed difference is probably even larger. Estimates based on our model solution to the Fokker-Planck equation (Section 2.3) suggest that the central power to the beam ions is twice as large with counter-phasing.

Predicted neutron enhancement
To study the effect of the RF power absorption on the volume averaged neutron rate, we use a simple model for the beam ion distribution function patterned after Stix's classic treatment of fundamental minority heating [33]. The model (Appendix) ignores spatial variations in the RF power deposition and beam distribution. Anisotropy in velocity space is approximated using formulas suggested by Hammett [34]. The Fokker-Planck equation for the distribution function f is solved in terms of the energy W ; the most important terms in the equation are the beam source, the collisional drag on electrons and thermal ions, and the angle averaged quasi-linear RF diffusion operator [35], where K RF is a constant that depends on the wave field and is proportional to the RF power and D RF is proportional to the Bessel function J 2 l−1 (k ⊥ ρ i ). (Here inj Figure 11. Model beam ion distribution function versus energy for three different values of k ⊥ ρinj. The parameters are the ones used to fit the data in Fig. 10. k ⊥ is the perpendicular wavenumber and ρ i (W ) is the fast ion gyro-radius.) This simple model reproduces the observed parametric dependences of the neutron rate. Within the framework of the model, the RF power deposition K RF is a free parameter. (The beam deposition rate also is a free parameter but does not affect the predicted neutron enhancement.) All the other parameters in the model are measured. The RF parameter K RF is adjusted to yield the average neutron enhancement for the parameters of discharge 96043. Then, with all other parameters held fixed, the dependence of the distribution function f on k ⊥ ρ inj is computed. The results are shown in Figs 10 and 11. The predicted neutron enhancement agrees well with the measured values (Fig. 10). The calculated distribution function both below and above the injection energy is quite sensitive to the value of k ⊥ ρ inj (Fig. 11). For k ⊥ ρ inj = 1.1, the distribution function is essentially a slowing down distribution and no enhancement to the neutron emission is predicted or observed. Near k ⊥ ρ inj = 1.5, a tail begins to appear and the neutron rate is enhanced. The enhanced neutron emission is generated by tail ions above W inj and by the 'hump' of additional ions just below W inj , with roughly equal contributions from each region. For fourth harmonic heating at this Alfvén speed, the maximum value of RF diffusion coefficient D RF occurs when W 250 keV (Fig. 12), so any increase in ρ inj results in larger RF acceleration.
The model also successfully predicts weak central fifth harmonic heating. For the present experimental parameters, the RF diffusion coefficient at the injection energy is 2.2 times smaller at the fifth harmonic than at the fourth harmonic (Fig. 12). Using the same parameters as in Fig. 10, the predicted neutron enhancement for central fifth harmonic heating is only 1.12, in excellent agreement with the measured value of 1.14 ± 0.05.

Effect on stability
In this section, enhanced sawtooth stability (Section 3.1), Alfvén instability (Section 3.2) and the eventual sawtooth crash (Section 3.3) are discussed.

Enhanced sawtooth stability
Under some conditions, combined beam and 60 MHz RF heating doubles or triples the sawtooth period ( Fig. 13(a)). If enhanced sawtooth stability occurs, it virtually always appears immediately after onset of the RF (in the first or second sawtooth cycle). When the plasma conditions remain approximately constant, many enhanced periods are observed, as in Fig. 13. In cases where only modest increases in period are observed (as in Fig. 2), the first or second period is usually longest; presumably, the modest (º25%) increase in density associated with the RF increases the collisional drag on the beam ion tail, decreasing the energetic ion density below the threshold for sawtooth stabilization.
In many cases, the large 'monster' sawtooth crash after a period of enhanced stability causes a transition to H mode (presumably because the large heat pulse triggers a bifurcation in edge transport). An example is shown in Fig. 14. Immediately after the RF is turned on, the sawtooth period doubles from ∼100 ms to ∼200 ms. The neutron rate is ∼50% larger than the classical prediction during this time. The first 'monster' crash at 1.97 s triggers an H mode transition (evidenced by the characteristic reduction in D α light). Since the antenna is operated at a fixed peak voltage, the reduced antenna loading [17,36] causes an immediate reduction in coupled RF power. Meanwhile, the plasma density begins to increase in response to the improved particle confinement. Initially, the transition has little effect on either the beam ion tail (as indicated by the neutron enhancement) or the sawtooth period. In fact, the tail is sustained with ∼50% of the initial RF power! Eventually, however, on the timescale of the density rise, the neutron rate and sawtooth period approach the levels observed in the absence of the RF. Presumably, in the high density H mode plasma, collisional drag dominates over RF acceleration. Enhanced sawtooth periods were never observed in the 1998 campaign if the RF pulse began after the H mode transition.
Enhanced sawtooth periods are not observed in the absence of a beam ion tail; however, some plasmas with a large tail do not have enhanced sawtooth periods (Fig. 15). Apparently, tail formation is a necessary, but not sufficient, condition for enhanced sawtooth stability. The large scatter in Fig. 15 is not  caused by differing positions of the resonant layer, as most of the plasmas with large values of neutron enhancement were centrally heated. In L mode, approximately 1 MW of coupled 60 MHz power is needed to affect the sawtooth period appreciably (Fig. 16). If sawtooth stabilization was already established during the L mode phase, lengthened sawtooth periods sometimes persisted at lower For the H mode data (squares), the coupled power was typically a factor of two larger in the L mode phase (crosses) of the discharge (cf. Fig. 14). power levels (0.5-0.7 MW) during the H mode phase (the H mode points with T ST > 200 ms in Fig. 16).
Enhanced sawtooth stability occurs when the fourth harmonic resonance layer is near the magnetic axis (Fig. 17). To date, long sawtooth periods have been observed only for beam powers of 2.5-3.3 MW (not for >4 MW or <2 MW). Stabilization is not observed forn e ² 4 × 10 13 cm −3 . The sawtooth period does not correlate with plasma shaping parameters such as the elongation or the gap between the outer wall and the plasma (correlation coefficient |r| < 0.2). The dependence on the edge safety factor has not been studied yet; all discharges have q 95 =4.6-5.2.

Alfvén modes
Alfvén activity is often observed in plasmas with lengthened sawtooth periods. Frequently, a mode is first detected by the magnetic diagnostics when the central electron temperature saturates (Fig. 18); the mode then persists until the sawtooth crash.
These persistent Alfvén modes are probably TAEs. Figure 19 shows the gap structure of ideal MHD [37] for the discharge shown in Fig. 18. The relatively strong magnetic shear causes significant radial variation in the position of the toroidicity induced gap. This indicates that a mode of large global extent will suffer continuum damping [38,39] or mode conversion [40]. Because the beam injection is unidirectional, the DIII-D plasma rotates in the direction of the plasma current. This implies that the mode frequency in the plasma frame ω pl differs from the measured frequency in the laboratory frame ω lab by the Doppler shift, ω pl ω lab − nω rot , where n is the toroidal mode number. Since the radial location of the Alfvén mode cannot be determined from the magnetics measurements, the toroidal rotation profile ω rot is used to estimate ω pl as a function of radial position [41]. The measured mode frequency falls in the TAE gap near the centre of the plasma and in the gap induced by ellipticity [42] (the EAE gap) near the plasma edge. Since the Alfvén activity correlates strongly with central quantities such as T e (0), the V2 neutral particle signal, and the neutron rate (discussed later), it is reasonable to assume that the mode is excited in the plasma interior. If this assumption is correct, the persistent Alfvén modes are TAEs.
The gradual evolution of the mode frequency in Fig. 18 is probably caused by changes in the background plasma parameters. The frequency drops 13%  Figure 19. Continuum of ideal MHD versus normalized radius ρ at 2.02 s in discharge 96043. The measured frequency including a local correction for the Doppler shift (solid line) is shown; a toroidal mode number n = 5 is assumed. The q profile (dashed curve) is derived from the EFIT reconstruction that fits magnetics and MSE data. The dotted curve indicates the higher frequency mode shown in Fig. 21. between 1.915 and 2.025 s. During this time, the electron density at ρ = 0.2 increases 25%, which should lower the TAE frequency ∼12%. The q profile also evolves ( Fig. 13(b)), which shifts the radial location of a given mode. In addition, beta increases 16%. Theoretically, pressure effects tend to raise the frequency by elevating the entire gap structure [37] but tend to decrease the frequency by shifting the mode from the centre of the TAE gap towards the lower continuum [43]. We conclude that the measured changes in n e and β can account for the gradual evolution of the measured frequency and are consistent with identification as a TAE. [On the other hand, the number of fast ions (as inferred from the neutrons) also increases 15%, so conceivably energetic particles also affect the frequency.] The persistent TAEs only occur in discharges with long sawtooth periods (Fig. 20) that have more than ∼1 MW of coupled 60 MHz power. Usually a single peak is observed, although occasionally a doublet is observed. The toroidal mode number ranges from n = 4 to n = 7, with n = 6 typically observed. In unstable beam heated discharges in DIII-D, the TAE toroidal mode number is usually n = 3-5 [44] while, during RF heating in larger tokamaks [45][46][47], n 10 is typical.
In some discharges, higher frequency Alfvén activity is observed (Fig. 21). In the plasma frame, these modes have the frequency of an EAE in the plasma interior and of a triangularity induced Alfvén eigenmode (NAE) at the plasma edge (Fig. 19). The higher frequency modes are found in plasmas with smaller changes in the sawtooth period (Fig. 20); another difference is that they often appear after a sawtooth crash and disappear later in the cycle  Figure 20. Laboratory frequency of Alfvén activity (from magnetics) versus duration of the sawtooth period for all of the 60 MHz data in the 1998 campaign. Discharges without any detectable activity (crosses) are assigned zero frequency. Some discharges have modes that persist until the sawtooth crash (squares) while, in other discharges, the modes usually appear after the sawtooth crash (diamonds). In H mode, the spectrum is complicated (asterisks) by ELMs and by modes that appear even in the absence of the RF power.  Fig. 2. The high frequency mode that appears after the sawtooth crashes is probably aliased at the Nyquist frequency of 500 kHz. The contour intervals are the same as those in Fig. 18.

AE FREQUENCY (kHz)
E A E TA E V A /4πqR (kHz) Figure 22. Frequency of Alfvén activity versus the nominal TAE frequency for the persistent (squares) and high frequency intermittent (diamonds) data in the 1998 campaign. The measured frequency is corrected for the Doppler shift using ω pl = ω lab −nωrot, where ω lab is measured by the magnetics, n = 5 is assumed and ωrot is the central toroidal rotation measured by the CER diagnostic. The data in the box are corrected for aliasing (about the Nyquist frequency of 500 kHz). The abscissa is calculated from the vacuum toroidal field and the line averaged electron density; q = 1 is assumed. The lines indicate the expected central TAE and EAE frequencies, respectively.
( Fig. 21). The frequency of these modes gradually increases during the sawtooth cycle (Fig. 21), but this is probably caused by aliasing (see below); most likely the true frequency gradually decreases as the density increases, as for the persistent TAEs. The data are from all of the discharges on one runday (Figs 2-5). The V2 line of sight was excited by a Right source beam (crosses) and by a Left source beam (squares); the axes are shifted (but not rescaled) to lie on the same graph. The diamagnetic signal is normalized to agree with the equilibrium signal following the RF heating (as in Fig. 3). Discharges without detectable Alfvén activity are assigned zero frequency. The horizontal error bars represent typical uncertainties.
Another difference between the persistent TAEs and the intermittent higher frequency modes is that the central electron temperature saturates in discharges with TAEs ( Fig. 18) but does not saturate in discharges with higher frequency modes (Fig. 21).
The density and toroidal field dependence of the mode frequency are consistent with identification as Alfvén eigenmodes (Fig. 22). The central frequency of the persistent modes agrees well with the expected frequency of TAEs in the plasma core (correlation coefficient r = 0.86). For most of the discharges, the Nyquist frequency was 500 kHz. The frequencies of the intermittent higher frequency modes apparently are aliased. After reflecting the data about the Nyquist frequency, the data agree well with the expected EAE frequency (correlation coefficient r = 0.80). On a few discharges, data were acquired at 2 MHz but no EAE or NAE activity in the range 500-1000 kHz was observed.
Neutral particle data are available for the discharges with EAE activity. The modes are observed only in the plasmas that have a large efflux of high energy neutrals from the centre of the plasma (Fig. 23(a)); Alfvén activity is absent in discharges with little increase in central flux during the RF pulse. Also, the discharges with EAE activity are more anisotropic than the stable discharges ( Fig. 23(b)). Evidently, a centrally located, anisotropic, energetic beam ion tail is required to excite Alfvén modes. However, the EAEs often appear immediately after a sawtooth crash (Fig. 21), so they may be excited when the centrally heated energetic ions are redistributed to the plasma periphery.
The onset of persistent TAE activity apparently causes T e to saturate. The correlation noted in Fig. 18 is found in all of the discharges with persistent TAE activity (Fig. 24). Saturation of T e (0) is only observed in plasmas with Alfvén activity. Even though the central electron temperature saturates, the volume averaged neutron rate continues to increase until the sawtooth crash (Fig. 18). In fact, it is likely that the Alfvén activity begins to affect T e (0) even before the modes are detected by the magnetics; on TFTR, internal Alfvén modes are observed by a reflectometer diagnostic ∼20 ms before they are detected by edge magnetic probes [47]. Presumably, the Alfvén modes cause anomalous beam ion transport in the plasma core. The redistributed beam ions continue to produce neutrons, but the central beam ion density saturates near the marginal stability point [48]. With the central beam ion density fixed, the power transfer from the energetic ions to the electrons saturates, which results in saturation of T e (0). Changes in electron transport associated with the TAE activity may also play a role.

Sawtooth crash
Eventually sawtooth stability is lost and a 'monster' crash occurs. An example of the sawtooth crash is shown in Fig. 25. There is a 6 kHz n = 1 magnetic precursor that grows at a rate of ∼3 ms −1 . The crash itself appears as a spike on several magnetic channels; the relative delay of the spike varies rapidly with toroidal angle. This particular crash triggers a persistent n = 2 successor oscillation. The bandwidth of the radiometer ECE signal was limited by the electronics to ∼25 kHz; the crash duration is shorter than 0.02 ms. The island size measured by ECE is very small ( 1 cm) up to the last period prior to the crash. In contrast, during beam heating alone, the island grows to ∼15 cm prior to the crash. A rapid 15% reduction in neutron emission occurs at the sawtooth crash.
The MSE measurements indicate that, with RF heating, the crash occurs for a lower value of q 0 and a larger value of the q = 1 radius than with beam heating alone (Figs 13(b,c)). Detailed analysis of the q profile just before and after the crash at 2030 ms is shown in Fig. 26(a). Within experimental uncertainty, q 0 returns to unity after the sawtooth crash, indicating that complete reconnection of poloidal flux inside the q = 1 radius must be taking place. The MSE pitch angle measurements (converted to an equivalent poloidal field profile) are shown in Fig. 26(b), where the change in slope from before to after the crash is clearly visible. Note that there are nine MSE channels inside the q = 1 surface, which gives high confidence that the change in q 0 during the sawtooth cycle is accurately determined. The uncertainty in q 0 due to statistical fluctuations in the MSE data is small, of the order of δq = ±0.01, while systematic errors are estimated to be ∼±0.05. The effect of the radial electric field [25] on the MSE measurements has been considered and is negligible in these discharges. The inferred q = 1 surface is in good agreement (within 2 cm) with the inversion radius measured by ECE. Relaxation of q 0 1 was reported previously for DIII-D [49], but these changes in safety factor (∆q 0.2) are much larger in this case. In many instances, partial reconnection events are observed within 50 ms of the monster sawtooth crash (cf. T e traces in Figs 2,8 and 13). Although the gross features of the q profile are consistent with a flat q profile inside the inversion radius immediately after the sawtooth crash, these subsequent partial sawteeth suggest the existence of some non-monotonic fine-scale structure in q(r).
The normalized q = 1 radius is always approximately 0.30 at the time of the crash during the RF heating ( Fig. 13(c)). This constancy of the q = 1 radius is observed for all of the L mode discharges with TAEs in the 1998 campaign (Fig. 27), even though the duration of the sawtooth period varies by a factor of two! Irrespective of the rate of central current diffusion, once the q = 1 radius reaches 0.30 ± 0.02, a crash occurs. Alternatively, the data can be interpreted as a limit on the central safety factor to values above 0.83 ± 0.02.
The q = 1 radius at the time of the crash is ∼50% larger than that with beams alone (Fig. 13), presumably because the accelerated tail ions provide additional sawtooth stability. Theoretically, the number of energetic ions N b needed to stabilize the sawtooth depends sensitively on the radius of the q = 1 surface r 1 , scaling as r 3 1 [12]. Because the TAE activity causes saturation of the central beam-ion population, discharges with similar Alfvén gap structures have similar values of N b , even though the RF absorption varies. Since the TAE activity enforces the condition N b const, the discharges all lose sawtooth stability at the same value of r 1 .
In H mode, the limiting value of r 1 is smaller (∼0.20), but the limiting value is still independent of the sawtooth period (Fig. 27). The difference between L mode and H mode plasmas could be caused by changes in the q profile that alter the ideal MHD stability of the sawtooth (e.g., the magnetic shear near the q = 1 surface may change) or it may be caused by changes in the Alfvén gap structure that lower the threshold for TAE activity.

Discussion
Fourth harmonic heating of beam ions is an effective heating scheme in these experiments. The incremental energy confinement time is usually larger than the confinement time prior to the RF heating (during beam heating alone), presumably because the fast ions are better confined than the thermal plasma.
The density was modest in these experiments. If the coupled power is increased to overcome increased collisional drag (P RF ∝ n e at constant temperature), fourth harmonic heating should also work in higher density plasmas. Because of the dependence of k ⊥ on the Alfvén speed, higher density increases k ⊥ ρ i (∝ √ n e ), which leads to stronger beam ion absorption (Fig. 10). On the other hand, in plasmas with high ion temperatures, second harmonic heating of thermal hydrogen could absorb most of the power. The dependence on antenna phasing (Fig. 8) is similar to observations on JET [50]. Since differences in direct electron heating efficiency are also observed in DIII-D [51], it seems likely that the stronger beam ion heating with counter-current propagation is caused by differences in wave propagation [52] rather than by a particle pinch [50]. The shift of the resonance layer relative to the magnetic axis caused by the Doppler effect may also play a role.  [11], TFTR [66] and JT-60 [5]. (Adapted from Kimura et al. [5].) The enhanced sawtooth stability seems consistent with the conventional picture of precessional drift stabilization [12]. Deeply trapped ions with energies of 80 keV precess at a frequency of ω d cW/(ZerRB T ) 1.6 × 10 5 rad/s, which is an order of magnitude higher than the observed Doppler corrected sawtooth precursor frequency. For moderate anisotropy, the poloidal beta of the fast ions β ph must exceed πα 1 where α 1 is a constant that depends on the details of the distribution function, 1 is the aspect ratio of the q = 1 surface, β MHD p is the poloidal beta in the absence of fastion stabilization and β p is the achievable poloidal beta with fast ions [12]. Using the volume averaged poloidal beta before and during RF power injection to estimate β MHD p and β p , respectively, the expected threshold in β ph /α 1 is ∼0.08, which is of the same magnitude as the observed anisotropy of the plasma (Fig. 3).
In terms of density and safety factor, the operational space for sawtooth stabilization is similar to that observed with fundamental minority and second harmonic heating schemes in other large tokamaks (Fig. 28). The required RF power is smaller, but this is primarily due to the smaller volume of DIII-D. The threshold power density of ∼50 kW/m 3 with fourth harmonic beam ion heating is comparable to the threshold power density for second harmonic beam ion heating in JT-60.
The TAEs are probably excited through resonances with the precessional motion of the trapped beam ions. Theoretically [53], the resonance condition is ω T AE = nω d + mω b , where ω d and ω b are the bounce and precessional frequencies, respectively, and n and m are integers. If the mode is driven by deeply trapped 100 keV ions, the required toroidal mode number is n ω T AE /(ω d + ω rot ) 7, which is consistent with the measured value.
Unlike the complicated 'tornado' spectra in JT-60U [54], the 'chirping' spectra in TFTR [47] and the 'pitchfork' spectra in JET [46], the Alfvén spectra in our experiments usually consist of a single steady mode. (Occasionally two modes briefly coexist.) Differences in power probably account for the different behaviour: in our experiment, the RF power is barely adequate to stabilize the sawtooth or destabilize the TAEs, while the JT-60U, TFTR and JET plasmas are driven well above threshold. Gradual saturation near the marginal stability point is readily achieved in a weakly driven system [55], but more complicated phenomena occur in strongly driven systems. (It should be noted, however, that, although there is a general tendency in the TFTR data for more peaks to appear in the spectrum with increasing RF power [48], complicated spectra are sometimes observed near the threshold for sawtooth stabilization [47,56].) At the sawtooth crash, our absolutely calibrated MSE diagnostic indicates that q 0 returns to unity (Fig. 13). The temporal behaviour of q 0 is similar to data obtained on JT-60U during second harmonic heating [57]; these data were corroborated by the TAE behaviour. Complete reconnection is consistent with the Kadomtsev model [58] of the sawtooth crash.
The large 15% reduction in neutron emission at the sawtooth crash (Fig. 25) is not surprising. The central density drops roughly 10%, which causes some reduction in the beam-plasma rate, but the beam ions are probably also redistributed. Theoretically, redistribution is expected when the crash duration is shorter than the precessional period [59]; in our case, the crash duration (<20 µs) is significantly shorter than the precessional period of full energy beam ions (∼40 µs).
What is required to stabilize sawteeth indefinitely? Naively, one would increase the RF power to increase the number of tail ions N b inside the q = 1 surface, but this is unlikely to work. Both sawtooth stability and TAE instability are caused by trapped beam ions. The saturation of the electron temperature (Fig. 12), together with the steady, gradually growing amplitude of the TAE, strongly suggest that the central beam ion density saturates at a level that is determined by the stability properties of the TAE. The situation is illustrated in Fig. 29.
Driving the system harder with additional RF power only results in stronger oscillations about the TAE threshold beta, but does not enhance sawtooth stability. In fact, oscillations about the marginal stability point trigger the sawtooth crash at a slightly lower value of β p (MHD drive) than would otherwise be achieved, as suggested by the observation that 'chirping modes' trigger monster sawteeth in TFTR [47]. The key to enhanced sawtooth stability is to raise the TAE threshold. The TAE threshold can be increased either by increasing the TAE damping or by decreasing the TAE drive. A sheared gap structure can increase TAE damping [44,60,61], but the gap structure (Fig. 19) is already relatively narrow in these plasmas. Theoretically [43,62], increasing the plasma beta can also increase the damping. Alternatively, one can modify the fast ion distribution function to maximize sawtooth stabilization while minimizing TAE drive. Theoretically, for both the internal kink and the sawtooth, the fast ion term in the dispersion relation is of the form [63] γ where the diamagnetic frequency ω * h is proportional to the gradient of the hot particle pressure and the brackets represent an average over the distribution function f . For the TAEs, the mode frequency ω is typically much larger than the precessional frequency ω d but, for the sawtooth, ω ω d . This implies that the ratio of destabilizing TAE drive to stabilizing sawtooth drive scales as According to Eq. (4), the TAEs weight the high energy portion of the distribution function more heavily than the internal kink. This implies that the ideal distribution has many particles at intermediate energies to provide sawtooth stabilization, but little 'tail' at high energies to destabilize the TAE. To achieve the ideal distribution function, one chooses the injection energy W inj to correspond to the lowest possible precession frequency that provides sawtooth stabilization. For high harmonic heating, the RF diffusion coefficient D RF has a maximum value for a particular value of k ⊥ ρ inj (Fig. 12). The harmonic number and Alfvén speed are selected so that the maximum of D RF occurs near the injection energy. This gives maximal acceleration near W inj , which produces a 'hump' in the distribution function. Above W inj , D RF decreases (the acceleration is weaker), so tail formation is minimized. The RF power and density are adjusted so K RF τ se is modest; if this parameter is too small, RF acceleration is negligible while, if this parameter is too large, an excessively long tail is generated. (Of course, in a reactor, other constraints such as the alpha particle birth energy may dictate the effective injection energy.) Even with N b optimized through TAE control, a sawtooth crash will eventually occur if the q = 1 radius grows too large (since N b ∝ r 3 1 ). Unless there is some form of current profile control, the same thermal instability that drives the sawtooth cycle will eventually overwhelm the additional stability margin provided by the fast ions. With non-inductive current drive to tailor the current profile, sawtoothfree operation with q 0 º 0.8 should be possible. For DIII-D, off-axis electron cyclotron current drive in the direction of the plasma current [64] is an attractive possibility.

Conclusions
Fourth harmonic heating of deuterium beam ions creates an energetic perpendicular beam population inside the q = 1 surface. The absorption increases as k ⊥ ρ inj increases. Second harmonic absorption on thermal hydrogen is important for concentrations above ∼5%.
The precessing beam ions can transiently stabilize the sawtooth when the RF power exceeds 1 MW and the resonance layer is located inside the q = 1 surface.
The trapped ions also destabilize the TAEs and EAEs. The TAEs are responsible for saturation of the central electron temperature.
Sawtooth crashes are caused by current accumulation on-axis. Because the beam population is limited by the TAEs, the crashes always occur when the q = 1 surface reaches the same radius.
Complete reconnection at the sawtooth crashes occurs.
In future work, we plan to measure the beam ion distribution above 100 keV. These data are needed for quantitative comparisons with theories of wave damping, sawtooth stabilization and TAE destabilization. The dependence of sawtooth stabilization and magnetic reconnection on plasma shape will be investigated. We will also attempt to stabilize sawteeth for the duration of the RF pulse in a highl i plasma with high confinement [65] through noninductive current profile control.

Model distribution function
The beam ion distribution function f is a function of position r and velocity v. Since the neutron data are volume averaged, we seek a local solution and ignore the spatial dependence of f . The energy W and pitch angle ξ = cos −1 (v /v) are selected as velocity space co-ordinates. Since the beam-target neutron rate only depends upon energy, our goal is to express the distribution function in terms of the energy alone, f = f (W ).
The Fokker-Planck equation is where C represents the collision operator, Q represents RF acceleration and S represents sources and sinks. We seek a steady state solution, ∂f /∂t = 0. In