Particle distribution modification by low amplitude modes

Modiﬁcation of a high energy particle distribution by a spectrum of low amplitude modes is investigated using a guiding center code. Only through resonance are modes effective in modifying the distribution. Diagnostics are used to illustrate the mode–particle interaction and to ﬁnd which effects are relevant in producing signiﬁcant resonance, including kinetic Poincar´e plots and plots showing those orbits with time averaged mode–particle energy transfer. Effects of pitch angle scattering and drag are studied, as well as plasma rotation and time dependence of the equilibrium and mode frequencies. A speciﬁc example of changes observed in a DIII-D deuterium beam distribution in the presence of low amplitude experimentally validated ToroidalAlfv`eneigenmodesandreversedshearAlfv`eneigenmodesisexamined indetail.Comparisonwithexperimentaldatashowsthatmultiplelowamplitude modescanaccountforsigniﬁcantmodiﬁcationofhighenergybeamparticle distributions.Itisfoundthatthereisastochasticthresholdforbeamproﬁle modiﬁcation,andthattheexperimentalamplitudesareonlyslightlyabovethis threshold.


Introduction
Energetic ion populations often drive Alfvén waves unstable in toroidal magnetic confinement devices [1].Alpha particles or high energy beam ions may drive similar instabilities in ITER [2].A goal of current research is to predict and control fast-ion transport in ITER and other future devices.In particular, because ripple increases rapidly near the field coils, ripple loss of alpha particles and the resulting thermal wall load is a sensitive function of the assumed alpha particle density profile [3], so modification of the birth profile must be accurately predicted.To that end, quantitative understanding of fast-ion transport in existing devices is essential.
Recent experiments in the DIII-D tokamak are well suited for tests of theory.In these plasmas, deuterium beam ions drive many toroidicity-induced Alfvén eigenmodes (TAEs) and reversed shear Alfvén eigenmodes (RSAEs) unstable [4].The mode structure is measured with electron cyclotron emission (ECE) and beam-emission spectroscopy (BES) diagnostics.Careful comparisons of the ECE measurements with the linear eigenfunctions calculated by the MHD code NOVA show excellent agreement in the mode shape [4] and temporal evolution [5,6].Saturated mode amplitudes are derived by scaling the prediction of a synthetic ECE diagnostic applied to NOVA calculated eigenfunctions.Likewise, the same scaling factor gives quantitative agreement with the electron density fluctuations measured by BES [4], confirming that the mode amplitudes are accurately determined.The resultant beam-ion transport is measured by five independent techniques [7,8], including spatially resolved fastion D-alpha (FIDA) spectroscopy.The data imply strong central flattening of the fast-ion profile during the early phase of the discharge when many Alfvén modes are unstable [7,8].In plasmas without appreciable MHD activity, the FIDA profiles agree well [10] with the profiles predicted by the NUBEAM module [9] of the TRANSP code but, in the presence of the strong TAE and RSAE activity, the profile in the inner half of the plasma is much flatter than classically expected [7,8].
Since the Alfvén activity and fast-ion transport are correlated [7,8], it is natural to assume that they are causally related.To test the assumption that the Alfvén modes cause the additional fast-ion transport, in previous works [7,8], we inserted the magnetic part of the NOVA calculated eigenfunctions that were experimentally validated by ECE measurements into the guiding center code ORBIT [11] and calculated the expected fast-ion transport.A similar approach was successfully applied to observations of fast-ion transport by fishbones in the poloidal divertor experiment [12] and tearing modes in numerous toroidal devices [13][14][15][16].There have been many previous numerical studies of fast-ion transport by Alfvén eigenmodes [17][18][19][20][21][22][23] but, in all cases, appreciable fast-ion transport occurred when the mode amplitude was δB/B ∼ 10 −3 .However, for the DIII-D Alfvén eigenmode case, this procedure failed, underestimating the observed transport by an order of magnitude.Apparently, in contrast to the transport by a few large-amplitude (δB/B ∼ 10 −3 ) low-frequency modes, the same approach is inadequate for many small-amplitude (δB/B ∼ 10 −4 ) modes in the Alfvén frequency band.The goal of this work is to understand why previous attempts to simulate DIII-D failed [7,8].
We find that the failure was due to neglect of the electrostatic part of the perturbation, important only for modes with frequencies well above those of typical fishbones.In general, many small amplitude Alfvén eigenmodes can cause fast-ion transport that approaches the experimentally observed levels, and simulations can reproduce this provided that all modes and all important effects are included in the simulation; that is, the guiding center equations must include many harmonics and all significant mode-particle coupling terms.The beam-ion transport possesses a stochastic threshold very near the experimental mode amplitude values, so the results are very sensitive to small effects.Even including the electrostatic potential, truncating the spectrum by omitting the smallest harmonics results in failure to reproduce the results.
To be able to analyze different cases, we investigate specific harmonics to find their ability to modify the particle distribution through the production of small scale islands in the phase space of the particle trajectories.Mode amplitudes are taken to be the experimental value of the order of δB/B 2 × 10 −4 , and at this amplitude modes that are not resonant have no effect on the distribution.Small scale islands contribute to particle transport in the presence of pitch angle scattering by producing excursions from the nominal drift surfaces and there is the possibility of stochastic transport if there are a sufficient number of islands large enough to cause overlap.Other possible contributions to distribution modification are slow changes in resonance location, caused by the slowing down of the high energy particles on the electrons, mode frequency sweeping and equilibrium modification.In section 2 we review the guiding center equations in general magnetic coordinates including the effect of a flute-like magnetohydrodynamic mode.Next we examine many possible mechanisms for the rapid modification of a high energy particle distribution.Section 3 examines the mode-particle resonances capable of changing a high energy particle distribution, and the diagnostic tools useful to find the effect of particular modes.Section 4 derives a means of finding the effect of a time dependent equilibrium on the distribution, section 5 the effect of the potential produced by electrons eliminating the electric field parallel to the magnetic field and section 6 the additional complication of high beta modes with compressional components.Section 7 treats the effects of pitch angle scattering and drag.Next, in section 8, we examine the case of a particular discharge in DIII-D, and examine these processes to determine those that are effective in causing fast-ion transport.We find that the system possesses a stochastic threshold, and that the experimental values are only slightly above it, making the simulation very sensitive to mode amplitudes and other small effects.Section 9 presents the conclusions.

Guiding center equations
Using the contravariant and covariant representations for the equilibrium field B, Here 2πψ p is the poloidal flux and 2πψ is the toroidal flux between the axis and the surface given by ψ p = constant, dψ/dψ p = q(ψ p ) and 2πI is the toroidal current inside ψ, 2πg is the poloidal current outside ψ and δ is a measure of the nonorthogonality of the coordinate system.The coordinates ζ and θ are toroidal and poloidal angles, respectively.We use straight field line coordinates, in which the local field helicity is a function only of the flux surface.It has been shown that the projection of the particle orbit in the poloidal plane and the toroidal precession are independent of δ [24].We will thus set it to zero in equations having to do with particle orbits.This modification also makes the Lagrangian equations of motion Hamiltonian in character.Figure 1 shows a sample DIII-D equilibrium, obtained by solving the equilibrium code ESC [25] (based on the experimental equilibrium found using EFIT [26]).
The equilibrium evolves during a 30 ms period and the minimum q value evolves in time from 4.312 to 4.145.Data from the equilibrium code are then put in spline form and used in the guiding center code ORBIT [11].Unstable modes in the plasma are represented by time dependent perturbations of the equilibrium field.The simplest perturbation of the equilibrium is a flute-like perturbation, with δ B = ∇ × α(ψ p , θ, ζ, t) B. The equations of motion are easily modified to account for the effect of α [11,27], which has the dimensions of length, and has a Fourier representation α = m,n α mn (ψ p ) sin(nζ − mθ − ω n t). ( This perturbation produces magnetic islands at resonant surfaces where q(ψ p ) = m/n.The equations of motion are given by the Hamiltonian form where the Hamiltonian is H = ρ 2 B 2 /2 + µB + , is the electric potential, and the canonical coordinates and momenta are θ , ζ and P θ , P ζ , with Here ρ = v /B, the normalized particle velocity parallel to the magnetic field.Distance is normalized to the radius of the magnetic axis and time is normalized to the on-axis cyclotron frequency.Plasma rotation is easily included by making a function of the flux surface.
The guiding center equations including flute modes are [27] where Here we use the notation f α to denote ∂ α f .The code ORBIT [11] consists of a fourth order Runge-Kutta implementation of these equations along with routines for particle insertion, collisions and diagnostics.

Resonance
Only through resonance with a perturbation is a significant modification of the particle distribution possible, since the mode amplitudes are known to be very small, with δB/B 2 × 10 −4 .Particle resonance is not equivalent to magnetic field resonance.There is significant drift modification from a simple field line following analysis, in which ζ = q θ , and the perturbation is fixed in time.Toroidal precession considerably modifies ζ at high energy.
Particle resonances are determined not by the q profile of the magnetic field, but by an analogous 'kinetic q factor' determined by the particle motion, including the effects of precession and the large shift of the drift surfaces [28].For low particle energy, zero frequency and pitch λ = v /v = 1 the kinetic q factor approaches the magnetic q profile, but precession and shift of the drift surface strongly depend on frequency, energy and pitch.We are interested in passing particle resonance.The θ and ζ time dependence can be written in the form [27] where ω t = v /R is the transit frequency and ζd , θd are drift terms second order in ρ/R with ρ the cyclotron radius.
A qualitative understanding of resonance can be gained by examining the large aspect ratio circular equilibrium case.Using B = 1 − r cos θ gives the large aspect ratio expressions ζd (ρ B 2 +µ)r cos θ, θd −(ρ B 2 +µ) cos θ/r, the θd term dominates the energy evolution and one finds [27] for the particle energy evolution in the presence of a single harmonic with Q m = nζ − mθ − ω n t.Thus the harmonic m produces two surfaces at the points where Q m−1 and Q m+1 are resonant due to the cos(θ ) dependence of the drift terms.
It is fairly easy to assess the effect of a particular mode on the particle distribution by examining a Poincaré plot for a particular choice of either co-moving or counter-moving particles, which we refer to as a kinetic Poincaré plot to distinguish it from a plot of the magnetic field.Points are plotted in the poloidal cross section whenever nζ − ω n t = 2πk with k integer.Neglecting the effect of the drift in modifying the toroidal motion then gives successive Poincaré points, with ζ = ω t t satisfying For there to be a periodic fixed point in θ with period m we also require θ = 2πl/m with l integer.But we also have θ = ω t t/q, giving this last equation determining the location of the resonance.Note that the poloidal mode number m does not appear in this expression.A resonance appears whenever there exist integers m , l such that this relation can be satisfied.Thus for q > m /ln resonance occurs with a co-moving passing particle and for q < m /ln it occurs for a counter-moving passing particle.Since the Alfvèn frequency is generally large, only rapidly moving particles are capable of participating, and important interaction occurs only for high energy heating particles or for fusion products such as alpha particles.Note that for co-moving passing particles (ω t > 0) and n > 0, increasing the mode frequency ω n increases the q value of the resonance, and increasing energy or pitch (and thus ω t ) decreases the q value of the resonance.These islands exist in real space and in the energy variable.The excursions in energy and flux due to the mode are related through [27] However, equation (11) was calculated neglecting the drift motion and for a large aspect ratio circular equilibrium.In order to examine the effect of resonance for arbitrary energy and which can be solved by Newton's method.The deposition energy of particles in the plot is fixed at the magnetic axis and then determined at other flux surfaces by equation (13).It typically varies by 10% or 20% over the minor radius, decreasing outward from the magnetic axis.
A kinetic Poincaré plot shows islands indicating resonance of the particles with the perturbation, and it includes all nonlinear couplings, deviation of the orbits from a single flux surface due to drift, and particle precession rates.The examination of kinetic Poincaré plots can easily be done in any equilibrium.
Figures 2 and 3 show kinetic Poincaré plots for different modes in an equilibrium with the reversed shear q profile of figure 1, to illustrate the islands produced by modes with δB r /B of order 10 −4 .Energies and pitch are chosen to reflect values near the peak of a DIII-D deuterium ion beam distribution.The dependence on energy, frequency and m values is shown.Note that the number of islands in the chain is not simply given by the mode number m, the value it would have by following field lines.Following the analysis above, the number of islands in the chain is denoted by m .All these plots have values of µB/E approximately equal to 0.6, a typical value for the DIII-D beam distribution.The large displacement in ψ p versus θ of approximate cos(θ ) form is the drift motion of the co-injected beam, found by solving for the unperturbed orbit using E = ρ 2 B 2 /2+µB with ρ = (P ζ +ψ p )/g with both P ζ and E constant.Figure 2 demonstrates the resonance position shift of an m = 10 island chain due to a change in energy, decreasing the energy increases the q value, and for this radius the surface moves outward.This resonance was produced by a single harmonic with m = 9 and n = 3.The first plot of figure 2 and the first plot of figure 3 produced by a single harmonic at 58 kHz but with different m values show that the resonance position and the value of m are independent of the m value.The two plots of figure 3 show the frequency dependence of the resonance surface, motion to a larger q value for increasing frequency.For this reversed shear profile the inner resonance moves inward and the outer resonance outward.Note that island size is proportional to the square root of the local magnitude of the perturbation, given by the eigenmode structure of the instability.Thus a shift in position of a resonance can also mean a substantial change in island width.
Another means of finding the ability of a particular mode to produce change in the particle distribution is to look for time averaged energy transfer to or from the particles [29], which can only happen if a particle is trapped in an island produced by the perturbation.Figure 4 shows plots of the P ζ , µ plane for an energy of 25 keV in the equilibrium shown in figure 1.
Here ψ w is the value of ψ p at the last closed flux surface.Each point in this plane corresponds to a particular orbit.The rightmost parabola symmetric about P ζ = 0 consists of orbits that pass through the magnetic axis.The small parabola centered about P ζ /ψ w = −1 consists of orbits that contact the last flux surface at the inner midplane, and the larger parabola consists of orbits that contact the last flux surface at the outer midplane.In the left plot particular orbit types are labeled.Domain A consists of confined counter-passing orbits.In domain B are the confined trapped orbits, and in domain C the confined co-passing orbits.Here co and counter refer to the direction of the plasma current.Domain D consists of counter-passing stagnation orbits, they do not circle the magnetic axis, and execute a small circle on the high field side of the axis.Domain P consists of so-called potato orbits, the parallel velocity is always positive, but they circle the magnetic axis by virtue of drift motion.Most other large domains are particle orbits that are promptly lost.For a derivation of these domains and further discussion see [27].
Initially a distribution consisting of all confined orbits is launched, and then a mode is allowed to act on the distribution.A time average is taken, and only particles showing a significant time averaged energy change are plotted.These plots show bands of resonance, indicating the particular orbits that are involved.The plots show resonance occurring for deeply and for barely passing particles in domain C, giving also the flux surface location.Whereas the Poincaré plots of islands show resonance only for a particular value of pitch, these plots indicate all pitch values giving resonance at one energy, appearing as a separate band for each value of µB 0 /E.Note that the 10/3 mode has a much smaller effect on the distribution than the 9/2 mode, there being no orbits which experience a 10% energy change in the latter case.

Time dependent equilibrium
For most of the simulations presented here, we use numerical equilibria obtained from discharge data, giving a good representation of the actual plasma equilibrium in the device.But each such realization is for a particular time in the discharge, and it is not presently feasible to incorporate time dependence in this manner.To investigate the role of a time dependent equilibrium we use instead a large aspect ratio low beta equilibrium, B = ∇φ + I (r)∇θ, I = r 2 /q, where the coordinate system is given by the minor radius coordinate r, the poloidal angle θ and the toroidal angle φ, with the major radius equal to X = 1 + r cos(θ ), q the field helicity dφ/dθ and the toroidal flux is given by 2πψ with ψ = r 2 /2.Here B is normalized to the value on the magnetic axis and the distances are normalized to the major axis.This analysis could easily be generalized to include an equilibrium at higher beta with a definite Shafranov shift.
Using this equilibrium, we can make the q profile time dependent.In the simplest form, in an axisymmetric configuration, but including a toroidal electric field and a time dependent q profile with the toroidal field constant but the poloidal field time dependent, the guiding center equations for the particle motion become where the denominator D = q + I + ρ I .Here B is the magnitude of the magnetic field and the electric potential.A prime indicates a derivative with respect to ψ p , the poloidal flux, given by The last term in the φ equation describes plasma rotation.
For any Hamiltonian variable the total time derivative is ṗk = −∂H /∂q k + ∂ t p k .The only variable thus modified through the ∂ t term by the changing equilibrium is ψ p .Without the ∂ t ψ p term the particle orbits erroneously convect inward or outward depending on the sign of ∂ t q.
Time dependence is introduced through q = q(ψ, t) with the toroidal field and hence toroidal flux constant in time, so we find directly the convective term This simply keeps ψ and thus r constant.
To simplify the numerics use a polynomial representation for the inverse q profile 1 q(ψ, t) A term linear in r is not possible, since this would make dq/dψ p infinite at the axis.Then since ψ p = dψ/q, we find analytically ψ p (ψ, t) and ∂ t ψ p (ψ, t).For the inverse, an initial guess of ψ = ψ p /i 0 followed by a Newton iteration works very well, since q is always positive and thus ψ and ψ p are monotone functions of each other.The primary effect of a time dependent equilibrium is the slow change produced in the location of mode-particle resonances.However, this effect is not different in kind from that produced by mode frequency modulation.In the DIII-D case studied this effect is very small and can be neglected.

Potential due to zero E .
Ion motion is further modified by the fact that the rapid mobility of the electrons makes the electric field experienced by the ions parallel to the magnetic field equal to zero.Thus in general it is necessary to add an electric potential to cancel the parallel electric field induced by d B/dt, with In Boozer coordinates, used in our simulations, taking = m,n m,n e i(nζ −mθ−ωt) the solution is but in general coordinates where I = I (ψ, θ) the solution is complicated by the coupling of different poloidal harmonics.The contributions to energy change and to motion across flux surfaces are proportional to (ρ/R)α from the magnetic perturbation and to (ω/ω 0 )qα m,n /(nq − m) from this potential, where ρ is the cyclotron radius, R the major radius and ω 0 the cyclotron frequency.Ideal modes vanish at the rational surface, so α m,n (ψ p ) is zero at nq = m.However, the potential is not zero at this point, and thus can have a significant effect near rational surfaces, provided ω/ω 0 is comparable to ρ/R. 4igure 5 shows kinetic Poincaré plots with and without the electric potential given by equation ( 22) for a 50 kHz TAE mode and 25 keV deuterium ions, giving in the DIII-D equilibrium ω/ω 0 = 3 × 10 −3 , ρ/R = 10 −2 .The effect of the potential can be neglected only for modes with ω/ω 0 ρ/R.If these terms are comparable the potential can have a significant effect on island size.This potential was not included in our original simulations [7,8].

Compressional modes
In a high beta equilibrium Alfvén waves take on a compressional component.The modification of the guiding center equations upon the introduction of a flute-like mode such as used in the previous sections, with δ A ∼ B is well known.Motivated by the coupling of the lowfrequency shear Alfvén spectrum with the acoustic spectrum [34,35] and the fact that both the RSAE [33] and other modes can have a significant compressional component, we now wish to find the modification of these equations upon introducing a mode with a general compressional component, with δ A ⊥ B. Write the Lagrangian in the form with ρ = v /B the normalized parallel velocity, H = ρ 2 B 2 /2 + µB + the Hamiltonian and the electric potential.Introduce a perturbation δ A = a θ ∇θ + a ζ ∇ζ and require B • A = 0, giving δ A = a(ψ p , θ, ζ )[q∇θ − ∇ζ ] and δ B = (aq) ψ p ∇ψ p × ∇θ − a ψ p ∇ψ p × ∇ζ + (qa ζ + a θ )∇ζ × ∇θ .Make the usual Fourier decomposition a(ψ p , θ, ζ, t) = a mn (ψ p ) sin(nζ − mθ − ωt).Also δ A is time dependent, so it introduces an electric field with giving Lagrange's equations are d dt and noting that da/dt = a ψ p ψp + a θ θ + a ζ ζ + ∂ t a we then find the equations of motion with denominator and f = I a ψ p +g(aq) ψ p , where C = 1−ρ g ψ p +a ψ p , F = q +ρ I ψ p +(aq) ψ p and W = a θ +qa ζ .
The Hamiltonian is also modified by the introduction of δ A, since to first order ( B +δ B) 2 = B 2 + 2 B • δ B and thus we find δB = B 0 f/Z, with Z = gq + I , where we have used ∇ψ p • (∇θ × ∇ζ ) = B 2 0 /(gq + I ) and dropped terms of order a 2 , and B 0 is the equilibrium field.Thus Note that the canonical momenta are now P ζ = ∂ζ L = ρ g − ψ p − a, and P θ = ∂ θ L = ρ I + ψ + aq, and the terms in ∂ t a in equation ( 27) produce the electric field which is perpendicular to B. The magnitude of the compressional vector potential A is proportional to the plasma β = P /(2B 2 ).For a TAE mode we use [30] δB B −1 [∇P • ξ ⊥ + γ P ∇ • ξ ] with P the plasma pressure.This expression serves to determine f/Z and thus the function a mn (ψ p ), hence proportional to the plasma β.For the DIII-D discharge studied, with β = 1%, the terms in the guiding center equations are too small to be relevant, but could be important in high β discharges.

Collisions and slowing down
The collision operator for pitch angle scattering on a background is given by λ = λ(1 − νdt) ± (1 − λ 2 )νdt (30) with time step dt, collision frequency ν and λ = v /v.This is the numerical implementation of a Lorentz collision operator [31].Although it is typically very small, and by itself produces a small modification of the high energy profile through neoclassical processes, in conjunction with the phase-space islands produced by the perturbation it plays an important role in the simulations, allowing particles to diffuse between nearby resonance surfaces.As seen both in the kinetic Poincaré plots and in the plots of the P ζ , µ plane showing energy transfer, the resonance domains can occupy a very small domain in µ or in ψ p , and thus a small amount of pitch angle scattering or change in the velocity magnitude can cause a particle to move in or out of resonance.
The slowing down due to electron drag is simply given by where the drag frequency ν s must be calculated for the density and energy of ions and electrons under consideration.As shown in figure 2, decreasing the energy increases the q value of the resonance.Thus the major effect of the slowing down of the distribution is a slow movement of the resonance surfaces.Particles trapped in resonance with the mode can thus be carried inward or outward for some distance until they escape due to scattering or a loss of resonance.But there can be an even stronger motion of these surfaces caused by changes in the mode frequency, shown in figure 3, with motion to a larger q value for increasing frequency.Similarly, if the equilibrium is changing on a slow time scale, the modification of the q profile causes motion of the resonance surfaces.The effects of slowing down, equilibrium change and frequency chirping can cancel or add, depending on the sign of the frequency change.These effects should be included together if they are known.They can all produce a form of bucket transport [32] through the motion of the islands.

General considerations
A quantitative comparison between theory and experiment is challenging.One source of difficulty is that the fast-ion diagnostics measure quantities that depend on the fully evolved fast-ion distribution function f .In the experiment, f develops on the slowing-down timescale τ s , which is much longer than the timescale of wave-particle interactions with the Alfvén waves.On this longer timescale, beam fueling and Coulomb collisions are important.The competition between these processes in shaping f can be written as an evolution equation, where S is the source term, C represents Coulomb collisions and W the effect of the waveparticle interactions.In the experiment, all three terms on the rhs are of comparable importance in the determination of f , some with more importance in some parts of phase space than others.For quantitative comparison with experiment, the computed f is employed in forward modeling of the expected signals that include reaction cross sections, spatial and velocity-space resolution, etc (see, for example, [10]).
Quantitative comparison with theory is relatively simple if the instability causes a sudden, transient change in the distribution function, for example, for a strong bursting instability such as a fishbone.In this case, the wave-particle interaction dominates the evolution of f for a short time t τ s .One can use a code such as TRANSP to compute f before prior to the burst, use a code such as ORBIT to analyze the wave-particle interaction during the burst and calculate the subsequent distribution function f after , then use forward modeling to compute the expected change in measured signals M/M before .Unfortunately, the situation is far more complicated for quasi-steady transport, such as the Alfvén eigenmode case considered here.Presumably, in the actual experiment, the many Alfvén modes sweep through the distribution function as they are driven unstable, causing transport where the phase-space gradients are largest.There is no sudden change in f that can be detected by the temporal evolution of a diagnostic signal.In a complete theoretical treatment the real distribution function f would drive the modes unstable, the nonlinear dynamics would be computed and f would relax due to properly modeled transport in phase space.Meanwhile the source would replenish f so that it drives other modes unstable, with the entire complicated process followed for a full slowing-down time.Unfortunately, a realistic treatment of f over such a long timescale is too expensive for current computer codes.Quantitative comparison with experiment is difficult.
Out of expediency, we employ a simpler approach.The NUBEAM module in TRANSP properly computes the source and collision terms S and C but does not treat the wave-particle interactions W .To estimate the phase-space averaged effect of the wave-particle interactions, we have employed an ad hoc diffusion coefficient D B in NUBEAM.Through trial and error, the magnitude and spatial profile of D B is adjusted to yield a distribution function that is consistent with the experimental measurements [8].This provides a quantitative estimate of the fast-ion transport but the details in phase space are probably incorrect.To improve this, with ORBIT, we select the Alfvén modes at a particular instant in time.We use harmonics given by NOVA.We run the simulation for times that are short compared with τ s and use the change in the distribution function f to estimate a phase-space averaged change in f .With this procedure, we find that the predicted transport is comparable to the level observed in the experiment.The details of this comparison for shot 122117 at t = 340-370 ms [6] follow.
Note that for high energy particles the flux surface is not a good descriptor of an orbit.As seen in figure 2 the drift motion is a significant fraction of the minor radius.TRANSP is axisymmetric, so the particle distribution in it can be written as f (E, P ζ , µ).During the action of the modes, E and P ζ are not constant, on a short time scale they either oscillate about the initial value or are trapped in a resonance.In addition, pitch angle scattering modifies µ.A guiding center simulation does a short time scale average over dt = of these variables for each particle.Besides the oscillations there is a longer time scale secular motion.It is necessary to construct the new f by making a grid in E, P ζ , µ and adding up particles in each domain.This gives the new f = f n .Then assuming linearity in time, we divide f n − f by , giving df/dt (E, P ζ , µ) for use in TRANSP.This df/dt cannot be characterized as flow or diffusion, it probably has elements of each.And as we will see, the motion could even be subdiffusive and nonlocal.The phase-space islands produce local flattening of the distribution on a time scale given by the trapping time in the islands and this flattening is not diffusive.This time can be easily obtained by looking at kinetic Poincaré plots at short time intervals, observing how many transits it takes for the points to circle the island O-points.Additionally, as the islands form and particles circle in them, they also scatter in and out of resonance, giving a diffusion process with a step size given by the island width.Also, if there is island overlap stochastic orbits could result in nonlocal transport.The combination of these effects uniquely define f n , giving a prescription for advancing the TRANSP distribution.
Note that to reconstruct the distribution given the values of E, P ζ , µ for TRANSP it is necessary to place all particles on the outboard midplane.This is the only place that all confined particles are guaranteed to pass 5 .Attempting to place particles back at their original poloidal angle after the change in E and P ζ would often result in no solution.

Beam particles
The initial beam profile was obtained from a TRANSP calculation of the distribution function [9], with energy ranging from 20 to 80 keV.TRANSP produces a list of 10 5 to 10 6 particles characterizing the deuterium beam, giving the velocity and position of all particles at a particular time.There are beam particles present with lower energy than this, but we are interested in the effect of the modes on high energy particles only, so the distribution is truncated at 20 keV.The energy distribution is shown in figure 6, along with the distributions in poloidal flux and in pitch, with the pitch expressed in terms of the magnetic moment µ and energy E. The pitch is λ = v /v = ± √ 1 − µB/E, with v the velocity.The distribution is almost entirely co-passing, and significantly peaked around µB 0 /E = 0.6, with B 0 the on-axis field strength.The distribution in energy has a dominant contribution at E 25 keV.The flux distribution gives the number of particles in equal size zones of the square root of poloidal flux; it is approximately proportional to the particle density.

Mode spectrum
A large spectrum of TAE and RSAE modes was observed to be present with amplitudes in the range of δB/B ≈ 10 −4 , as determined by density and temperature fluctuation measurements.The spectrum of modes we use in the simulation is given by NOVA, with the amplitudes of the various modes fixed by comparison with temperature fluctuation measurements.An example of the comparison of a NOVA calculated eigenmode with ECE data is shown in figure 7(a) for a f = 78 kHz TAE, where the perturbed electron temperature (δT e ) is plotted versus the normalized square root of toroidal flux (ρ) [6].For comparison with ECE measurements, a synthetic diagnostic as described in [4] was used to process the NOVA predicted temperature perturbation.The actual poloidal harmonic content/structure comprising the TAE is shown in figure 7(b), where it is seen that at least 10 harmonics contribute significantly, something typical of the global TAEs discussed here.By scaling the NOVA prediction using a single constant to match the ECE data (figure 7(a)) the amplitude of the perturbation wavefields is obtained.The inferred amplitude is shown in figure 7(c) as a function of major radius (R) along the device midplane, where the radial magnetic field perturbation (δB r ) is scaled to the local magnetic field strength B. For the majority of experiments on DIII-D, typically AE amplitudes obtained in this manner are found to be δB r /B < 10 −3 .The frequency dependence of the spectrum of modes included in the simulation is shown in figure 8.Only the RSAE modes have significant frequency variation over the range of time considered; the TAE modes are fairly constant in frequency.
To estimate magnitudes of the field components it is sufficient to use large aspect ratio approximations.The perturbed field is approximately giving components and it is the radial (across the equilibrium field) component that is most effective in producing modifications [27].Figure 9 shows the radial profiles of the relevant component of the field that directed across the equilibrium flux surfaces, for some of the harmonics, plotted along the outboard midplane.We used 133 harmonics in the simulation, more than 10 harmonics per mode, and a reduction of this number by dropping those with small amplitude significantly  One defect of the present analysis, which uses the eigenvalues for the instabilities obtained with the code NOVA, is that in general these eigenfunctions should change shape and location during changes in the q profile or during frequency chirping.Simulations of time dependent equilibria and frequency chirping with fixed eigenfunctions probably underestimate transport.In the actual situation the modes would track the changes in the distribution as well as equilibrium modification, always shifting to the location of maximum drive, and hence maximum profile modification.

Distribution modification
Figure 10 shows the result of application of the spectrum of figure 9, multiplied by two to give mean amplitudes of dB r /B of 2×10 −4 , to the beam distribution.The distribution modification is shown at 9.6 ms.There is a steady increase in the beam modification in time due to the slow effect of the interplay between the resonances and the pitch angle scattering.We also show the modification produced by pitch angle scattering alone, and by the waves alone.It is clear that there is a synergistic interplay between collisions and small phase-space islands.
For all values of collisionality and mode amplitude considered, the beam distribution modification is similar in shape to that shown in figure 10.To compare the effect of perturbations with different collisionalities and amplitudes, we introduce the mean distribution shift, through Scaling of this number with collisions and mode amplitude at a time of 4.8 ms is shown in figure 11.There is some induced transport for zero collisionality, with the value almost doubling from this ν = 0 value at physical values for this experiment (I = 0.008).The variation with perturbation amplitude is stronger, almost quadratic in dB/B.Recall that the resonance island size scales as the square root of the perturbation, so if this were the only effect the transport would scale as the island width squared, and hence linear in dB.The fact that it is stronger than this indicates that higher order islands due to nonlinear couplings are filling in the phase space, giving more rapid transport than what would be produced by simple resonance width increase.This explains why it is important to keep all relevant harmonics in the simulation.
Also shown in figure 11 are the results of time dependence of the RSAE modes, allowed to chirp during the simulation period.Two simulations were performed, at 10 and at 20 ms.Frequency chirping of the RSAE modes can produce a small increase or decrease in the profile modification, but appears not to be an important influence on distribution modification in this case.In this case, since upward frequency chirping moves resonances to larger q values, bucket transport actually hinders the observed change in the beam profile.
For a quantitative comparison, we note that δn/n(0) from figure 10 is 3% at ψ p = 0.3, corresponding to a radius of 0.3 m.Because of the steep density profile, δn/n locally is more than 10%.This is the correct location for the observed maximum density increase [8], and using the 9.6 ms time of this result we calculate an effective diffusion rate of D 0.1 × (0.3) 2 /(0.01) 7 m 2 s −1 , in rough agreement with the ad hoc diffusion coefficient inferred from the fast-ion measurements of ∼5 m 2 s −1 .Neoclassical diffusion alone, as seen from figure 10, gives about 0.02 m 2 s −1 , in agreement with usual diffusion estimates, giving about 0.01 m 2 s −1 .We conclude that the observed small-amplitude TAEs and RSAEs can account for the flattening of the fast-ion profile.This analysis does not yet lead to a full predictive theory of the effect of MHD modes on high energy particles, because we have had to use the experimentally determined mode amplitudes.To attain a complete analysis the nonlinear mode saturation must also be predicted.But what we have demonstrated is that low amplitude modes can be sufficient to produce large effects, and that at high frequencies any analysis must include the electric potential and the full spectrum of harmonics.
Figure 12 shows the time dependence of the beam modification.It is seen to be approximately linear in time out to 10 ms.This important result indicates that the process may be approximated using a fairly short time evaluation with the guiding center code, and extrapolated in TRANSP up to times when the distribution has changed sufficiently to warrant a new evaluation of the modes that are present.One time scale determining the process is the phase mixing time in an island.Numerically, we find that particles circle once in islands of this size in about 1 ms.Thus it takes several milliseconds for phase mixing in the islands to occur, and this could explain the apparent saturation of the modification at 10 ms.After phase mixing is completed the distribution can only evolve through island width dominated diffusion through collisions and the effect of stochastic orbits.Recall that we are performing an initial value simulation.In the actual experiment the continual resupply of ions by the source and the constant shifting of the modes due to distribution modification would eliminate this saturation effect.The rate of change in the distribution in this 9.6 ms simulation translates into a very significant modification of the beam profile during the 100 ms beam slowing-down time.
Equilibrium modification, at least in this experiment, is not a significant cause of beam profile modification.The change in the q profile during the course of the experiment produces only small changes in resonance position.

Stochastic web
It is not clear from the previous simulations whether the particles have large scale excursions in the plasma even in the absence of collisions, i.e. whether the mode perturbations produce a stochastic web allowing large scale transport.It could be that the profile modification was entirely due to the profile flattening caused by the islands and the subsequent collisional scattering among them, with no global stochasticity.A plot of island width and location for a given particle energy and a number of modes, as shown in figure 13, indicates that there is some overlap of islands for some minor radii and for some frequency intervals, but it is not feasible to produce such a plot including all harmonics in the simulation, and in any case such a plot, of necessity done for one energy and one harmonic at a time to clearly show island width, does not include nonlinear couplings which produce additional islands of smaller amplitude.Note that significant island width occurs only in the inner part of the plasma (ψ p < 0.5) and that the island locations move inward with increasing frequency.The plot shows the outboard midplane island width and location, so the mean position is even closer to the magnetic axis than shown.
To examine the question of chaotic transport we launch a distribution of 2000 particles all on the same drift surface, initially all with a single value of energy and pitch at the outer midplane, but distributed randomly toroidally.We choose a distribution characteristic of the beam, with a pitch of λ = 0.6 and an energy of 25 keV.Figure 14 shows the final particle positions after 7 ms in the presence of the modes used in the simulation, with magnitude dB/B 1.2 × 10 −4 (a) and dB/B 1.6 × 10 −4 (b) but no collisions.It is clear from this figure that the density level of the small scale islands in the simulation allows for global stochastic transport, with a stochastic threshold 1.2 × 10 −4 < dB/B < 1.6 × 10 −4 .In the interval between these numbers there is a gradual increase in the number of particles found outside the original drift surface after 7 ms, but it is clear that even at dB/B = 1.4 × 10 −4 the last KAM surface has been broken, and that a very long simulation would produce a flat distribution.Even for larger amplitudes the particle excursions are limited radially, indicating that the stochastic domain exists only in the core of the device, explaining the flattening of the beam distribution in the core, but the absence of changes further out in the distribution.This is apparent in figure 9, where it is seen that the modes have significant amplitude only in the plasma core.At an amplitude of 1.2 × 10 −4 a few particles have lost sufficient energy to become trapped, allowing them to move inward radially, but no particles move outward.This explains the failure of previous attempts to describe this profile flattening.A small reduction in the efficacy of the modes, or in the number of harmonics, moves the system below stochastic threshold, and the resulting particle redistribution is due solely to the slow collisional transport between island chains.
Figure 15 shows the time dependence of P 2 ζ , for the same ensemble of 25 keV particles with pitch λ = 0.6, again clearly showing the existence of a stochastic threshold for dB/B 1.4 × 10 −4 , with rapidly increasing diffusion for perturbation levels above this.Since the phase space of particle trajectories is not far above stochastic threshold, random phase approximations for the transport cannot be expected to be correct.In a similar situation in the reversed field pinch, where chaos was present but not at a level well above threshold, the transport was found to be subdiffusive [36,37], dominated by a spectrum of Levy flights, and described by a nonlocal Montroll equation.In a future publication we will report the results of applying the techniques used in those publications to the present problem.

Conclusion
We have considered means by which low amplitude modes can produce changes in a high energy particle population.Depending on the types of modes present and the frequencies involved, and the nature of the equilibrium, several effects are possible.At high frequencies it is important to include the effect of the electric potential.We have demonstrated that many TAE and RSAE modes can significantly modify a beam particle distribution, even with amplitudes of the level of dB/B 2 × 10 −4 .The relevant factor is the presence of islands in the particle trajectory space, for particles at energies and pitches characteristic of the distribution.The islands modify the distribution in two ways.The islands produce local flattening of the distribution, and if two island chains are in close proximity, stochastic orbits result.Even without this occurrence islands provide particle excursions from the initial drift surfaces, and in the presence of collisions the islands produce additional diffusion.Time evolution of the mode frequencies and the equilibrium can assist the beam profile modification by causing the resonance surfaces to move throughout the plasma volume, interacting with particles of different energies and pitches as they do so.Mode chirping can induce distribution changes through bucket transport in some cases, and plasma rotation must be taken into account if the background plasma is involved in the mode generation and the rotation rates are significant compared with mode frequencies.In this study we are looking at modes destabilized by the beam particles and at the effect of these modes on the beam, so plasma rotation is irrelevant.Pitch angle scattering and drag assist the transport due to the small resonance islands produced by the modes.Compressional effects of the modes are not found to be relevant at the beta values present, but could be important in some discharges.
We find that the modification of the beam profile in DIII-D from that predicted by TRANSP can be explained by the effect of the spectrum of low amplitude modes observed to be present in the discharge.Previous simulations failed because of neglect of the electric potential, important at high frequencies.The transport possesses a stochastic threshold, so it is very sensitive to small changes in mode content and amplitude.The fact that a sufficient number of perturbations can lead to stochastic transport in a Hamiltonian system is of course not new, it has been known since the proof of the KAM theorem and demonstrations with many models, and has also explicitly been demonstrated for TAE modes [38,39].Even at the low amplitudes present in the experiment the phase space of the trajectories is found to be stochastic, allowing slow but large scale modification of the distribution.The effect of the modes can be described as a phase-space dependent contribution to df/dt, to be added to the diffusion produced by Coulomb scattering induced neoclassical transport.This term can be much larger than the neoclassical contribution for some parts of phase space.
Although this analysis used the experimentally measured mode amplitudes and thus does not yet provide a fully predictive theory, what it does clearly demonstrate is that low amplitude modes, with a full spectrum of harmonics and complete wave fields taken into account, are able to produce significant modification of high energy beam-ion distributions, consistent with the DIII-D experimental observations.

Figure 1 .
Figure 1.The reversed shear DIII-D equilibrium, showing the decrease in the q profile during the 30 ms of the discharge we have studied.

Figure 2 .
Figure 2. Kinetic Poincaré plots for mode m/n = 9/3, showing energy dependence of the m = 10 resonances.Shown are the resonances for a 9/3 perturbation on particles of 25 and 24.5 keV.

Figure 3 .
Figure 3. Kinetic Poincaré plots for mode m/n = 10/3, showing frequency dependence of m = 10 resonances.Note that the resonant surface is unchanged from the m = 9 mode for 58 kHz.

Figure 4 .
Figure 4.A plot of the P ζ , µ plane, showing only 25 keV particle orbits with energy loss or gain due to the mode, for a m/n = 9/2, mode at 90 kHz (left, showing 10% change) and a m/n = 10/3 mode at 90 kHz (right, showing only 0.5% change).Both the co-moving and counter-moving parts of the distribution are shown (equilibrium of figure 1.)

Figure 5 .
Figure 5. Kinetic Poincaré plots for mode m/n = 9/2, showing the effect of the potential on 23 keV beam particles for a 50 kHz mode.

Figure 6 .
Figure 6.Initial Beam distribution in energy, pitch and flux surface.

Figure 7 .
Figure 7. (a) Synthetic ECE diagnostic prediction (solid) using NOVA calculated f = 78 kHz global TAE overlaid with ECE measurements (diamonds).NOVA prediction scaled by single constant to match ECE data.(b) Poloidal harmonics comprising TAE from panel (a).(c) Calculated radial component of magnetic field fluctuation along device midplane (versus major radius) using amplitude obtained from comparison with ECE data [6].

Figure 8 .
Figure 8.Time dependence of the modes in the spectrum.

Figure 9 .
Figure 9. Harmonics showing the radial field component.

Figure 11 .
Figure 11.Beam distribution modification versus collisionality (a) and mode amplitude (b).In the collisionality scan the amplitude is shown for dB/B = 2 and 3 × 10 −4 , and in the amplitude scan the collisions are fixed at ν = 2 s −1 ; the nominal values for the experiment are also shown with a larger point in each plot.The individual triangles and squares are results of RSAE frequency sweeping, performed at 10 ms (triangle) and at 20 ms (square).

Figure 13 .
Figure13.Island width and location, at the outer midplane, for a number of harmonics for beam particles at 23 keV and outboard pitch of 0.6 as a function of mode frequency.Island overlap is present at some frequencies for some radial locations.

Figure 14 .
Figure 14.Particle distribution after 7 ms, all launched at the outer midplane with pitch λ = 0.6 with random toroidal angle in the presence of the modes of section 8.3, but with no collisions.Amplitudes of dB/B 1.2 × 10 −4 (a) and dB/B 1.6 × 10 −4 (b).

Figure 15 .
Figure 15.Time dependence of dP 2 ζ , for different values of δB/B.