A description of the full-particle-orbit-following SPIRAL code for simulating fast-ion experiments in tokamaks Academic research paper on "Physical sciences"

lopscience

¡opscience.iop.org

Home Search Collections Journals About Contact us My IOPscience

A description of the full-particle-orbit-following SPIRAL code for simulating fast-ion experiments in tokamaks

This article has been downloaded from IOPscience. Please scroll down to see the full text article. 2013 Plasma Phys. Control. Fusion 55 025013 (http://iopscience.iop.org/0741-3335/55/2/025013)

View the table of contents for this issue, or go to the journal homepage for more

Download details: IP Address: 93.180.53.211

The article was downloaded on 30/08/2013 at 19:44

Please note that terms and conditions apply.

Plasma Phys. Control. Fusion 55 (2013) 025013 (23pp) doi:10.1088/0741-3335/55/2/025013

A description of the full-particle-orbit-following SPIRAL code for simulating fast-ion experiments in tokamaks

G J Kramer1, R V Budny1, A Bortolon2, E D Fredrickson1, GYFu1, W W Heidbrink2, R Nazikian1, E Valeo1 and M A Van Zeeland3

1 Princeton Plasma Physics Laboratories, PO box 451, Princeton, NJ 08543, USA

2 University of California-Irvine, Irvine, CA, USA

3 General Atomics, PO Box 85608, San Diego, CA 92186, USA

Received 25 July 2012, in final form 29 November 2012

Published 21 January 2013

Online at stacks.iop.org/PPCF/55/025013

Abstract

The numerical methods used in the full particle-orbit following SPIRAL code are described and a number of physics studies performed with the code are presented to illustrate its capabilities. The SPIRAL code is a test-particle code and is a powerful numerical tool to interpret and plan fast-ion experiments in tokamaks. Gyro-orbit effects are important for fast ions in low-field machines such as NSTX and to a lesser extent in DIII-D. A number of physics studies are interlaced between the description of the code to illustrate its capabilities. Results on heat loads generated by a localized error-field on the DIII-D wall are compared with measurements. The enhanced Triton losses caused by the same localized error-field are calculated and compared with measured neutron signals. Magnetohydrodynamic (MHD) activity such as tearing modes and toroidicity-induced Alfven eigenmodes (TAEs) have a profound effect on the fast-ion content of tokamak plasmas and SPIRAL can calculate the effects of MHD activity on the confined and lost fast-ion population as illustrated for a burst of TAE activity in NSTX. The interaction between ion cyclotron range of frequency (ICRF) heating and fast ions depends solely on the gyro-motion of the fast ions and is captured exactly in the SPIRAL code. A calculation of ICRF absorption on beam ions in ITER is presented. The effects of high harmonic fast wave heating on the beam-ion slowing-down distribution in NSTX is also studied.

(Some figures may appear in colour only in the online journal)

1. Introduction

In toroidal fusion devices the motion of the particles is determined by the Lorentz equations. The particle motion consists of the fast gyro-motion with typical frequencies in the MHz range for ions, the drift motion with typical

frequencies in the range of up to a few hundred kHz for ions, and the precessional motion which is typically on the order of a few kHz for ions. The fast gyro-motion is often neglected and the orbits are calculated at the center of the gyro-orbit because the inclusion of the gyro-motion makes the calculations very expensive in CPU time [1-4]. In the guiding

center approximation the magnetic moment of the particles, which is an adiabatic invariant [5,6], is taken to be conserved while the magnetic and electric fields and other relevant plasma properties are calculated at the guiding center location. The assumption made in guiding center theory is that the gyro-orbit is small compared with the characteristic scale length of the fields, a condition that is usually met for thermal ions in present day large tokamaks. Fast ions, however, can have substantially larger Larmor radii that can become comparable to the characteristic scale lengths of the plasma profiles, turbulence and magnetohydrodynamic (MHD) structures such as Alfven eigen modes and therefore, the guiding center

0741-3335/13/025013+23$33.00

© 2013 IOP Publishing Ltd Printed in the UK & the USA

approximation is not valid anymore [7]. This is especially the case in spherical tokamaks such as NSTX with their low toroidal fields. Another consequence of averaging over the fast gyro-motion in the guiding center approximation is that the interaction between ion cyclotron waves (ICRF) and the ions is not captured in guiding center calculations although efforts have been made to approximate the ICRF-particle interaction in guiding center codes with kick operators [8,9] but only a full-orbit treatment of the particle motion can capture the interaction between particles and ICRF.

With the advance of high-speed parallel computers it has become possible to follow the full orbits of fast ions in toroidal geometry and use statistically relevant ensembles of test particles to study the interaction between waves and particles without the limitations imposed in guiding center theory. The motion of charged particles in electromagnetic fields is governed by the Lorentz equations:

— = — (v x B + E),

— = vt dt

with m the mass, q is the charge, r is the position and v is the velocity of the particle while B and E are the magnetic and the electrical fields at the particle position.

Apart from the electromagnetic forces the fast ions are also affected by the thermal plasma through which they move. At high energies the fast ions suffer drag from thermal electrons that slow them down without a significant change in direction while less frequent collisions with thermal ions are responsible for changes in the travel direction of the fast ions, known as pitch-angle scattering.

In this paper, a description of the SPIRAL code is given and is interlaced with results from physics studies to illustrate the code capabilities. The SPIRAL code is written in the C language [10] and it makes use of the Numerical Algorithms Group (NAG) library routine: nag_ode_ivp_adams_gen (d02cjc) [11] which is a highly accurate, variable order, variable step size modified Adams method for solving the differential equations. When particles are followed in static magnetic fields only for ~50 ms, their energies are usually conserved better than 10-5. Due to the nature of the problem, following a large number of independent test particles on their paths through the plasma, parallel processing of individual orbits is straightforward and we have used the MPI package [12] to handle the parallel processing. The execution time increases linearly with the number of test particles and scales inversely with the number of CPUs on which the code is running. Great care has been taken that (partial) results are stored on disk as soon as they they are available and can be monitored during the run. When a computer failure occurs during a run, the SPIRAL code can very easily be restarted from that point, avoiding loss of CPU time. In the output file a large number of variables is written for access with a post-processor code so that the results can be studied in detail. Moreover, realistic first walls for different machines can be selected in the code so that particle losses and heat loads to the first wall can be studied accurately. The

calculation of the intersection of particle orbits with the wall are handled by the used NAG routine which will stop when the particle intersects the wall.

In section 2 the implementation of the equilibrium magnetic fields is discussed. Equilibrium solvers usually give the relevant fields on a discrete spatial mesh and a robust interpolation routine is presented that obtains the fields between the mesh points under the condition that Maxwell's equations are satisfied. In this section single particle orbits in NSTX and DIII-D are studied and it is shown that gyro-motion effects are significant in NSTX while they are small in DIII-D.

Slowing-down and pitch-angle scattering are discussed in section 3 and 3D equivalents of the 2D equations for pitch angle and slowing-down as given in [13] are derived.

For an accurate comparison between simulations and experimental results it is important to use realistic fast-ion distribution functions that reflect the ones present in the experiments. It is (almost) impossible to measure the full spatial and velocity distribution of the fast ions. Therefore, the SPIRAL code has the capability to either read a given fast-ion distribution or generate it internally as discussed in section 4 and a study in NSTX is presented on the effects of a tearing mode on the beam-ion distribution.

Toroidal-symmetry breaking fields are ubiquitous in tokamaks. They arise from the finite number of toroidal field coils, from ferritic steel that is needed for strength, or they can be induced deliberately with external field coils for shaping the plasma edge and suppressing edge-localized modes (ELMs) while losses induced by toroidally localized RF antennas can create hot spots [14]. In section 5 the numerical techniques are discussed that were used to implement those error fields in the SPIRAL code in such a way that Maxwell's equations are satisfied accurately. Instead of using a Fourier harmonic decomposition of the ripple fields we have chosen to specify the ripple fields on a cylindrical mesh so that highly localized perturbations, which are poorly described by a finite number of Fourier harmonics, can be included as well as conventional toroidal ripples. In this section, two cases are presented. In the first case the effects of a highly localized magnetic perturbation from a scaled mock-up Test Blanket Module (TBM) for ITER in DIII-D on the neutral beam ions and on the 1 MeV tritium population are simulated and compared with observations. In the second study the effects of resonant magnetic field perturbations (RMPs) on the beam ion confinement in NSTX are simulated to investigate whether the effects are large enough to observe RMP-enhanced beam-ion losses and/or changes in the confined fast-ion population that can be measured with the Fast Ion Da (FIDA) diagnostic [15].

In section 6 the implementation of MHD modes and RF waves in the SPIRAL code is discussed. MHD activity can have large effects on the fast-ion population via resonant interactions. Resonances between fast ions and a reversed shear Alfven eigenmode (RSAE) in DIII-D are studied. In addition, the effects of short (<1ms) bursts of torodicity-induced Alfven eigenmode activity, which are often observed in NBI-heated discharges in NSTX, on the neutral-beam slowing-down distribution are investigated (section 6.1).

One of the major heating schemes proposed for ITER is ion cyclotron range of frequency (ICRF) heating where RF

waves interact with the gyro-motion of the ions. Because this interaction involves the gyromotion of the particles it cannot be calculated with particle-orbit following codes that use the guiding center approximation which neglects the gyro-motion. Because the full gyro-motion is retained in the SPIRAL code, the interaction between the ICRF waves and ions appears naturally in the simulations. Two studies in ICRF absorption by fast ions are discussed in section 6.3. The first one is a study for ITER in which an estimate is calculated for the parasitic absorption of ICRF on the deuterium beam-ion slowing-down distribution (section 6.3). As a second example the interaction between high harmonic fast waves (HHFWs) and neutral beam ions in NSTX are studied. Because of NSTX's low magnetic field a large number of harmonic ICRF resonances are present in the NSTX plasma and particles can interact with several of those along their drift orbits. A preliminary study is presented of HHFW effects on a beam-ion slowing-down distribution showing that a significant fraction of the energetic beam ions is lost due to the HHFW (section 6.3). Conclusions are drawn in section 7. Where possible, comparisons between full orbit and guiding center calculations will be made but an exhaustive study between the two calculations is beyond the scope of this paper.

2. Equilibrium magnetic and electrical fields

For an accurate calculation of charged particle orbits in plasmas, the magnetic and electric fields should be known accurately in the region where the particle moves. The equilibrium magnetic fields in tokamaks are usually calculated on a discrete mesh with equilibrium codes such as EFIT [16] and q-solver. SPIRAL can read the output from those codes, which consist of the poloidal flux function, (R, Z), and the diamagnetic correction to the toroidal magnetic field. The radial and vertical components of the poloidal magnetic field are obtained by taking the derivatives of with respect to R, and Z:

B7 = -

1 d%(R,7) R d7 ' 1 d%(R, 7)

%(R, 7) = ££ aijTi(R)Tj(7).

1.5 1.0 0.5 0.0 -0.5 -1.0 -1.5

NSTX pulse

0.0 0.5

major radius [m]

0 0.5 1.0 1.5 2.0 2.5 time [|i s]

Figure 1. (a) A full orbit (red) and its guiding center orbit (blue) for a 90 keV deuteron in NSTX. The particle was launched at the LFS mid-plane parallel to the magnetic field. (b) The parallel and (c) the perpendicular velocities along the drift orbit.

computer accuracy over the entire region where the Chebyshev expansion is valid.

We will now study two single particle orbits in the equilibrium fields of NSTX and DIII-D and find that full-orbit effects are significant in NSTX while in DIII-D the effects are minimal. The particle orbit in a typical equilibrium field for NSTX is shown in figure 1 where a 90 keV deuteron, which is the maximum NBI injection energy in NSTX, was launched at R = 1.4m and Z = 0.0m with a pitch, v\\/v, of one which means for NSTX that the particle velocity is anti-parallel to the local magnetic field. (The parallel component, vy, of the particle pitch, vy/v, is parallel to the local magnetic field whereby the sign of vy is defined to be positive when the particle moves in the direction of the plasma current which is opposite to the toroidal magnetic field direction in NSTX). Both the full orbit (in red) and guiding center orbit (in blue) are shown, where the latter is obtained from the full orbit by evaluating the instantaneous guiding center which is given as

Rgc ^p

mv x B

Taking accurate numerical derivatives from a function that is given on a finite mesh is difficult because it has to preserve the divergence-free condition as dictated by Maxwell's equations which translates in this case to dRBR(R, Z) + dZBZ(R, Z) = 0 because of the toroidal symmetry. In SPIRAL we have solved this problem by fitting a 2D Chebyshev polynomial expansion [17] to the given discrete points of the poloidal flux function:

Once a good fit is obtained (usually 15-30 radial and vertical terms are needed), derivatives can be obtained analytically by rearranging the expansion coefficients (see appendix A). It can be shown easily that the magnetic field obtained in this way is divergence free as required by Maxwell's equations up to the

where B is taken at Rp, the particle location. There are two things noteworthy about this orbit. First, the gyro-frequency is only 12 times higher than the drift frequency which violates the assumption made in guiding center theory that the gyro-frequency is much higher than the drift frequency. Second, it can be seen from figure 1 that the full orbit has excursions of up to 3.7 cm away from the guiding center orbit which can be understood as follows: at the launch point the Lorentz force is zero because particle velocity is parallel to the magnetic field (figure 1(b)). Only when a perpendicular velocity component develops (figure 1(c)), a magnetic force arises and pulls the particle back onto its poloidal orbit.

When this orbit is calculated with the guiding center equations [18] the projection of the orbit onto the poloidal plane is virtually the same as the instantaneous guiding center

1.0 1.5 2.0 2.5 major radius [m]

time [|a s]

1.5 2.0 2.5 major radius [m]

0.0 x [m]

Figure 2. (a) A full orbit (red) for a 80 keV deuteron in DIII-D. The particle was launched at the LFS mid plane parallel to the magnetic field. (b) The parallel and (c) the perpendicular velocities along the drift orbit.

(equation (4)) with deviations between the two orbits of less than 1 mm. Notably in the guiding center calculation the parallel velocity of this particle remains constant along the orbit because its magnetic moment was chosen to be zero at the start. As a consequence, a poloidal transit of the guiding center in the guiding center calculation is completed slightly faster than one transit in the full-orbit calculation. This difference is reflected in the poloidal transit frequency which is 388.2 kHz for the full-orbit calculation and 398.6 kHz in the guiding center approximation, a difference of more than 10 kHz (or 2.6%). This difference is not caused by numerical differences in the magnetic fields because the same parametrization of the magnetic fields was used in the guiding center and full-orbit calculations.

The orbit for an 80 keV deuteron in DIII-D is shown in figure 2 where the particle was launched at the LFS mid-plane at R = 2.12 m with its velocity anti-parallel to the local magnetic field line (pitch: v^/v = 1.0 because the plasma current is opposite to the toroidal magnetic field direction). In this case the gyro-frequency is more than two hundred times higher than than the drift frequency and therefore, guiding center theory is more accurate than in the NSTX case. This is also reflected in the difference between the full orbit and guiding center orbit where the maximum difference between the two orbits is 1.2 mm which is indistinguishable in figure 2(a). The perpendicular velocity that the particle develops in this case is a very small fraction (less than 4%) of the initial parallel velocity as can be seen in figures 2(b) and (c) and therefore, the poloidal transit frequency difference between the full-orbit calculation (66.378 kHz) and the guiding center calculation (66.430 kHz) is small, 52 Hz.

Key-particle-orbit frequencies, such as the poloidal transit frequency, and toroidal transit and precession frequencies, which are quintessential for wave-particle interaction studies, are determined in the SPIRAL code using a two-dimensional Fourier transform of the poloidal and toroidal projection of the particle orbit. For the poloidal frequencies the radial

Figure 3. A full orbit (red) and its guiding center orbit (blue) for an 80 keV deuteron in DIII-D for which the frequency spectra in figure 4 were calculated. (a) Poloidal and (b) top view. The particle was launched at the mid-plane at x = 2.12m, y = 0.0m with a pitch of 0.8 and followed for one poloidal transition.

frequency [MHz]

60-1 0 1 frequency [MHz]

Figure 4. spectra of the poloidal (a) and toroidal (c) projection of the particle orbit shown in figure 3. with a zoom around zero frequency in (b) and (d). The poloidal bounce frequency, fb, is indicated in (b) while in (c) the toroidal transit frequency, ft, is shown.

and vertical coordinates are treated as real and imaginary quantities, respectively, for the complex fast Fourier transform. The output is a spectrum with the bounce frequency and its harmonics centered close to zero frequency while the cyclotron frequency and its harmonics form broad structures at higher frequencies. The broadening is caused by the varying magnetic field along the particle orbit. In a similar way the toroidal frequencies are obtained from the projection of the particle orbit in the mid-plane.

Poloidal and toroidal spectra for the orbit depicted in figure 3 are shown in figure 4. The orbit is for an 80 keV deuteron launched at the mid-plane at x = 2.12 m, y = 0.0 m with a pitch of 0.8. The fundamental cyclotron frequency is clearly visible between ±12 and 18 MHz while the second and third harmonic cyclotron harmonics are visible at twice

Figure 5. (a) Particle energy as a function of time for a deuteron that was launched at 1 MeV and (b) the slowing-down time calculated from equation (12). When the particle reaches thermal velocity the slowing-down is switched off in the SPIRAL code.

and three times the fundamental frequency, respectively. The poloidal bounce frequency at 49.6 kHz and its harmonics are visible in figure 4(b) while the toroidal transit frequency at 191.6 kHz is visible in figure 4(d) together with its harmonics.

Static electrical fields in plasma equilibria arise from plasma rotation and pressure gradients:

E(R, Z) = — V$ = Vplasma X B + Vp/q

Cf = VddA(1 - m/ - "2 dv(v + Vc3)f v2

with v = |v|, k = vy/v the particle pitch, v the collision rate given by

ne Zeff e2qf ln A

4ne2mimf v3

[20,21] and vd is the pitch-angle scattering rate given by

Vd = V' where the critical velocity, vc is calculated as

me /2Kb Te\3/2 4 mi y me

with vplasma the bulk plasma toroidal velocity, p is the plasma pressure, q is the charge of the fast ion and B is the equilibrium magnetic field. These static electrical fields can be included in SPIRAL calculations as an electrical potential function $(^p).

3. Slowing-down and pitch-angle scattering

Before discussing perturbative fields that are included in SPIRAL we will first discuss slowing-down and pitch-angle scattering in 3D geometry and the particle distributions that can be included in the code.

On their journey through the plasma the fast ions lose energy to the background thermal plasma. At high fastion velocities the energy loss is mainly to the electrons without substantial scattering but when the fast ions approach thermal velocity the interaction with the thermal ions becomes important and the ion interaction gives rise to pitch-angle scattering.

Traditionally, slowing-down and pitch-angle scattering are written in guiding center formalism [19] as a collision operator, C, for a distribution of energetic particles, f:

with me, (mi, mf) the electron (plasma, fast-ion) mass, qf is the charge of the fast ion and Te is the electron temperature. The first term of the right-hand side of equation (6) describes the pitch-angle scattering while the second term describes the slowing down whereby the v3 term is due to fast-ion-electron interactions while the vc3 comes from fast-ion-ion interactions. From the collision operator (equation (6)) the slowing-down equation:

dv (v3 + v^)

can be obtained. This equation is valid in both the guiding center formalism with v = (vy, vx) and in full-orbit formalism with v = (vx,vy,vz) and it is valid for velocities above the thermal velocity of the bulk plasma. When no precaution is taken the fast ions can slow down to zero velocity in simulations where the particles are followed for several slowing-down times as can be seen from equation (10). In long simulations one can avoid this unphysical situation by marking the particle lost to the thermal bulk when the particle energy is close to the local plasma temperature or by switching off the slowing-down acceleration are as follows:

dt v = —v

(v3 + vc3)

where ne is the electron density and the Coulomb logarithm, ln A, is given by A = 12nnekD with kD the Debye radius

with erf (x) the error function and Ef = 1 miv2 is the fast-ion energy. When particle energy is much higher than the thermal energy the error function is one and the particle slows down as usual as can be seen in figure 5(a) while near the thermal velocity the error function goes to zero and the de-acceleration is stopped. In SPIRAL one can either stop following the

particle when it reaches thermal velocity or following it to the requested end-time with the slowing-down force turned off as describe above.

The slowing-down time quantifies the time it takes to thermalize the particle and one can specify the velocity slowing-down time, rv, or the energy slowing-down time, te which is half the velocity slowing-down time. The velocity slowing-down time can be calculated from equation (10) as

|dtv| v(v3 + и;?) erf

(if - 1

kn = ko(1 - 2vdAi) ± У(1 - ko)2vdAt,

In the high energy limit (v ^ vc) this becomes ts = 1/v but when the velocity is near the critical velocity the energy slowing-down time shortens as is shown in figure 5(b). In SPIRAL measured electron temperature and density profiles are used. Therefore, the scattering rates and critical velocities are dependent on the location of the particle in the plasma.

Pitch-angle scattering in guiding center theory has been derived in [13] for an ensemble of fast ions. The new pitch, kn, is calculated from the old one, ko, at discrete time steps, At, as

where the time step is chosen in such a way that vd At ^ 1 and the ± sign is determined randomly at each time step. This operator fulfills the condition that |kn | < 1 — vd At. In fact, this constraint is slightly too restrictive because it can be shown that there is a small region in pitch space close to |kn| = 1 that cannot be reached by the particle. This region is given by |kn | > 1 — vd At and can be removed by ignoring the factor of two in the first term of equation (13). Removing this factor does not influence the random-walk process due to the pitch-angle scattering because the random process is governed by the second term in equation (13).

In contrast to pitch-angle scattering in the guiding center approximation with two velocity components, the full-orbit treatment requires scattering in three dimensions. In this section, we derive the full 3D pitch-angle scattering from geometrical considerations and obtain the 3D equivalent of equation (13). We will also derive the 2D pitch-angle scattering equation (equation (13)) without the extra factor of two in the first term.

In the 3D equivalent of pitch-angle scattering as depicted in figure 6, the fast ion is scattered over a small angle, AO, into a velocity cone around its original velocity, vold. The location on the cone, indicated with 0, is in this case the random variable and averaging is performed by taking a large number of small scattering time steps along the particle orbit.

For the 3D pitch-angle scattering we can choose without the loss of generality the following coordinate system: (p0 = vyo, p1 = vXo, p2 = po x p1) with vyo and vxo unit vectors defined by the particle velocity parallel and perpendicular to the magnetic field before scattering. The parallel and perpendicular scatter angles are now given by: AOy = AO sin(0) and AOx = AO cos(0). From figure 6 it can be seen that the new parallel velocity after scattering becomes

Figure 6. Pitch-angle scattering in 3D parallel/perpendicular space.

while the new perpendicular velocities along the p1- and p2-axis become

v± n = v sin(0 + A0y) cos(A0±), vip„n = v sin(A0±),

which corresponds to a rotation of the old velocity vector over an angle pair (AOy, AOx) and therefore, the new pitch, |kn|, after scattering is naturally restricted to equal or less than one.

Using now the fact that |AO | ^ 1 and expanding the geometrical functions in equations (14) and (15) we obtain the following expressions for the three velocity components: vyn = v(cos(O) cos(AOy) — sin(O) sin(AOy)) cos(AOx)

« (vyo(1 — AOj/2) — vXp1oAOy)(1 — AOi/2), vx n = v(cos(O) sin(AOy) + sin(O) cos(AOy)) cos(AOx)

« (vyoAOy + vxp1o(1 — AO|2/2))(1 — AOX/2), vxp2n = v sin(AOx) « vAOx. (16)

These equations are implemented in the SPIRAL code where the angle, 0, is drawn from a uniform distribution between 0 and 2n.

We can recover the Monte Carlo collision operator (equation (13)) from the 3D results of equation (16) by noting that the expectation of the parallel scattering angle is < AO2 >0 = AO2/2 where < ■ >0 stands for averaging over 0 along the full circle while averaging over the perpendicular scattering angle gives < AOx >0 = 0. After substituting those values in equation (16) we obtain for the parallel component:

= v||o(1 - AO2/4) ± v±noA0/V2,

where the ± sign which arises from the square root of AO2/2 acts as the random variable in the scattering process. If we now identify AO2/4 with vdt we obtain

щn = v cos(0 + AOy) cos(AOx)

vyn = vyo (1 - Vdт) ± У v2o2VdT,

<D 0.3

.c o t 0.2

0 5 10 15

time [ms]

Figure 7. Spread of pitch-angle distribution as a function of orbit following time. The solid line is a square-root fit to the data.

which is the same result as in equation (13) except for the factor two in the U||o term.

The spread in pitch of an ensemble of particles that is initialized at the same location and with the same pitch angle should follow a square-root behavior with time. Both the guiding center and the 3D formulation follow that behavior and are in excellent agreement with each other as can be seen from figure 7 where the spread in pitch of an ensemble of 2500 particles, initialized with zero pitch was followed in time.

There is an important difference in pitch-angle scattering between the guiding center and full-orbit formalism that can affect the calculated orbits. In the guiding center formalism, the orbit of the guiding center is continuous during the scattering while the particle orbit has (small) discontinuous jumps because v± changes discontinuously when the particle scatters. This can be seen from the expression for the gyro-radius, p = v^/mc. Those discontinuities in the particle orbit are physically unrealistic. In full 3D geometry, the scattering takes place at the particle location and therefore, the orbit of the particle is continuous which is physically correct while the guiding center orbit has discontinuous jumps at each scatter.

4. Particle distributions

Single-particle-orbit studies are useful to gain insight in the details of the particle-plasma interaction but statistically significant ensembles of particles are needed for fast-ion transport studies. A number of fast-ion distributions can be generated by SPIRAL such as a fusion alpha birth profile and two-temperature Maxwellian distributions. Fast-ion distributions can also be read from an external file. This option is very useful to import the distribution functions as calculated with other codes such as the TRANSP/NUBEAM [22,23] package.

Ensembles of fusion-born alpha particles can be generated by SPIRAL from a given alpha-particle birth profile which is a flux function. The test particles are drawn from that distribution and placed uniformly toroidally and poloidally on the flux surface. The particle pitch is drawn from a uniform pitch distribution, the alpha particle energy is set to its birth energy of 3.5 MeV, and the gyro-phase is drawn from a uniform distribution between 0° and 360°.

Maxwellian fast-ion distributions of test particles can be generated from a fast-ion density profile in combination with a parallel and perpendicular temperature profile whereby the profiles are given as flux functions. The parallel temperature can differ from the perpendicular one. In a similar way as fusion-born alpha particle ensembles, Maxwellian ensembles are drawn from the given fast-ion density profile and placed uniformly on its flux surface. The particle energy and pitch are now obtained from the requested local parallel and perpendicular temperatures while the gyro-phase is drawn from a uniform distribution between 0° and 360°.

The most flexible way of specifying the fast-ion distribution is to read the distribution from an external file on which the particle birth locations and velocities are written and this option will be exploited in simulations reported in this paper.

4.1. Effects of a tearing mode on a beam-ion distribution

The SPIRAL code has been used to calculate accurate statistical representations of the beam-ion slowing-down distributions in an NSTX discharge with and without the presence of a magnetic fluctuations. This study was motivated by the observation of a low-frequency mode that caused a substantial changes in measured FIDA signals. The changes in the FIDA signals are indicative of changes of fast ion distribution function. SPIRAL was used to determine the fastion redistribution associated with modes of tearing or ideal kink nature. Here we report results from the tearing mode case, while the detailed analysis of the ideal kink effect is discussed in [24].

The TRANSP/NUBEAM package was used to generate the beam-ion birth profile for the SPIRAL simulations from the three beam sources that operated at 90 kV acceleration voltage. In figure 8 the distribution of the Larmor radii at the particle birth time for the three energy components of the beam is shown. Larmor radii of up to 19 cm are present at birth which should be compared with the characteristic radial size of the mode (see figure 9(a)) which is on the order of 2 cm on the low-field side and around 10 cm at the high-field side. Because the Larmor radii are larger or comparable to the characteristic mode structure the guiding center approximation is not valid here and a full-orbit treatment of the wave-particle interaction is needed. (Only at thermal energies well below 1 keV the deuterium gyro-radius becomes comparable or smaller then the characteristic mode structure at the low-field side). Moreover, the Larmor radii are also significant compared with the characteristic scale length of the equilibrium field (~0.5 m) which is outside the validity of guiding center theory.

o t CC

—I-1-1-1—

—I-1-1-1—

—I-1-1-I—

—I-1-1-1—

5 10 15

Larmor radius [cm]

Figure 8. Distribution of Larmor radii for 90 keV beam ions at birth for a typical NSTX discharge. The full-energy component (red) consists of 43.3% the injected particles while the half energy component (blue) consists of 39.5% and the third-energy component (green) contains 17.1% of the injected particles.

In the calculations 30 000 particles were taken from the 3D beam-ion birth profile as calculated by the TRANSP/NUBEAM package and distributed uniformly in time between 0 and 25 ms. All of the particles were then followed until they reached 25 ms which is more than twice as long as the average energy slowing-down time. Two runs were made. In both runs the toroidal ripple, the radial electrical field from the plasma rotation, slowing-down, and pitch-angle scattering were included but only in the second run the tearing mode was included along the lines explained in section 6. We have used an m/n = 4/1 mode as shown in figure 9(a) with an amplitude of Br/B = 8 x 10-3 based on observations, and a mode frequency of 10 kHz. When the mode was taken into account the beam-ion losses increase from 17% to 23%. The calculated slowing-down distribution of the confined particles for the case without the mode is shown in figure 9(b) while in figure 9(c) effects of the mode on the distribution are taken into account. In figure 9(d) the difference between the distribution with and without the mode is shown. From this figure it can be seen that the mode mainly affects the distribution below 35 keV with pitches between 0 and 0.5 which corresponds to mainly trapped ions, while the distribution at higher energies and a more parallel pitch is hardly affected by the tearing mode. The FIDA system of NSTX which registered clear changes in the fast ion distribution in the experiment is not sensitive to the part phase space where the calculated changes occur due to the tearing mode. When a kink mode was used as was reported in [24] changes in phase space were found that were compatible with the FIDA measurements.

5. Static perturbed magnetic fields

Toroidal ripple and error fields are common phenomena in tokamaks and it has been shown that those fields can contribute

significantly to fast-ion losses [25]. Harmonic expansions in the toroidal angle are frequently used in particle-orbit following codes to include those ripple fields [26]. For ripples generated by a finite number of toroidal field coils such an approach works well but in some cases a toroidally and poloidally localized perturbation is created which cannot be decomposed adequately with a manageable number of harmonics. Such highly localized perturbations are expected in front of the TBM in ITER because of the significant amount of ferritic steel that is used in those modules [27]. Two examples of toroidal ripple fields are shown in figure 10 where a conventional ripple induced by a finite number of toroidal field coils is shown for NSTX while a localized ripple in DIII-D that was deliberately induced by a set of coils to create the expected (scaled) fields in ITER induced by one pair of TBMs. The fields in DIII-D were used to study the impact of such localized perturbed fields on the plasma performance [28] and first wall heat loads [29].

In order to include both the harmonic and localized ripples accurately in the SPIRAL code, the radial, bR(R, y, Z), and vertical, bZ(R, y, Z), components of the ripple field are given at a number of toroidally equally spaced planes on an unstructured mesh in (R, Z). On each toroidal plane a 2D Chebyshev polynomial fit is made for the bR and bZ while the toroidal component of the ripple field, by(R, y, Z), is computed from the condition that the magnetic field is divergence free. For the interpolation in the toroidal direction we have used a quadratic polynomial around each given toroidal plane:

bix(R, y, Z) = b'x(R, Z) + uX(R, Z)y + vix(R, Z)y2 (19)

with x either R or Z, and i the index of the ith plane. We demand further that the radial and vertical fields and their derivatives are continuous half way between two adjacent planes. In order to calculate the toroidal field we expand the functions b\(R, Z), uix(R, Z) and v\(R, Z) of equation (19) into a finite sum of Chebyshev polynomials and from the condition that the magnetic field is divergence free we can solve the unknown Chebyshev coefficients of the functions uix and vxi . This involves differentiating the expansions of b'x, uix and v'x with respect to R and Z, summing the various components, integrating with respect to y and making use of the facts that the field is periodic in y and that the Chebyshev polynomials are orthogonal. All those operations are performed analytically using the Chebyshev expansion coefficients and therefore, the ripple field is guaranteed to be divergence free in the calculations down to the numerical precision of the computer. In most ripple calculations it is sufficient to include about 10-15 coefficients in the radial and vertical direction to describe the ripple fields accurately as is shown in figure 10 for a conventional toroidal ripple in NSTX and a localized perturbation in DIII-D created by the TBM mock-up coils. In this figure it can be seen that the reconstructed toroidal field perturbation in SPIRAL is in excellent agreement with the one that was given.

Three examples of SPIRAL calculations with static magnetic fields are given here. In the first example heat

r e n e

—i-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1—

no mode present

pitch: Vy/V

¡^ 60 g

r e n e

with tearing mode

—i-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1—

0.0 pitch: V||/V

0.5 1.0 1.5 major radius [m]

¡^ 60 g

r e n e

-i—i—i—i—i—i—i—i—i—|-no mode with mode

■ _c

0.0 pitch: V||/V

Figure 9. (a) Radial magnetic field fluctuation of a tearing mode that was used to calculate the slowing-down distribution (b) without and (c) with the mode from the TRANSP/NUBEAM deposition profiles indicated in pink. (d) The difference between the slowing-down distributions without and with the mode present.

loads on the TBM surface are calculated and compared with measured temperatures. Full-orbit effects in this case only become apparent near the wall at the point where they are lost [29]. In the second example simulated triton losses are compared with neutron signals that were measured during the TBM mock-up experiments in DIII-D [28]. Full-orbit simulations are needed in this case because the tritons are born in the plasma at 1 MeV and have gyro-radii of up to 0.2 m in the DIII-D magnetic fields. In the third simulation the effect of RMPs on the beam-ion confinement is studied for optimizing measurable signals in the planning phase of NSTX experiments. Because of the low magnetic fields in NSTX only full-orbit calculations can give accurate predictions for those experiments.

5.1. TBM-induced heat load simulations

Heat loads on the carbon tiles that protect the TBM mock-up coils in DIII-D were calculated in the presence of the highly localized TBM error fields as shown in figure 11. Tile temperature rises were then calculated from the simulated heat loads and compared favorably with measured tile temperature rises. Details of this experiment are given in [28,29]. The heat load was caused by 60-80 keV beam ions that are lost shortly after injection before the ions have time to slow down significantly. In order to calculate accurately the location where the particle impacts the wall, it is important to take into account the full-particle orbit because the Larmor radius, in this case up to about 4 cm at the point of impact, is comparable or larger than the characteristic TBM field gradient scale lengths

90 180

toroidal angle [

80 40 0.0

-80 80

-80 80

-40 -80

; 1 -j

90 180

toroidal angle [

Figure 10. Toroidal magnetic ripple field for NSTX (left) and the TBM mock-up in DIII-D (right) both at the mid-plane low-field side plasma edge. From the given radial (bR) and vertical (bZ) components the toroidal component (by) is reconstructed (black line) and compared with the given field (red crosses).

near the tile surfaces. The heat loads from the SPIRAL simulations reproduce the measured temperature increase at the back of the 2.5 cm thick protective carbon tiles in the TBM experiments as was reported in [29].

5.2. TBM-induced triton losses

When NBI is injected in DIII-D, a small fraction of the fast deuterium beam ions react with the thermal deuterium plasma:

D D |3He + n(2.56MeV)50%,

Dbeam + Dplasma ^ Jp + t(1.01 MeV)50%.

By measuring the 2.56 MeV neutron signal one can deduce the 1.01 MeV triton birth rate in the plasma because of the equal branching ratios of the two reactions. In DIII-D those 1.01 MeV tritons have a large Larmor radii (figure 12), up to 20 cm, they are marginally confined, and the orbits of the confined tritons are very sensitive to perturbations of the magnetic field.

The confined tritons slow down and react with the thermal deuterium to form 4He and a 14.1 MeV neutron. Therefore, the ratio between the 2.56 and 14.1 MeV neutrons can be used as a sensitive measure for the triton confinement. In the DIII-D TBM experiments it was found that this ratio decreased by 33 ± 5% [28] when the TBM error fields were engaged.

For a SPIRAL simulation of the triton burn-up measurements we have simulated a typical DIII-D H-mode discharge (pulse 140153) with 5.8 MW of NBI power injected in which the TBM fields were applied. The triton birth profile was obtained from a TRANSP run for this shot. Two SPIRAL runs were made, one without and one with the TBM fields present. In both cases 5000 test particles were followed

260 270 280

toroidal angle [°]

Figure 11. Heat loads due to beam-ion losses on the protective tiles (outlined in yellow) of the TBM mock-up experiment in DIII-D as calculated with the SPIRAL code.

for up to 100 ms. The energy slowing-down time of those 1 MeV tritons is 175 ms, so at the end of the simulation the confined tritons had reached an energies of 600 keV or lower. In principle, the tritons should have been followed to around 70 kev, the energy where the DT fusion cross section peaks but for reasons given below the triton population that is left after 100 ms is on well confined orbits.

Many tritons are born on unconfined orbits and 53.8% are lost within the first 0.1 ms of the simulations as can be

DIIID pulse: 140153

in § 60

1= o o

SPIRAL simulation of DIIID shot 140153 -i—|—i—i—i—|—i—i—i—|—i—i—i—|

no TBM

30% loss

with TBM

40 60 time [ms]

1.5 2.0

Figure 12. A 1.01 MeV triton orbit (red) and its instantaneous guiding center orbit (blue) in DIII-D as calculated with the SPIRAL code.

seen in figure 13. The large majority of the particles that are lost immediately are born on trapped orbits with large banana widths. In fact, after 0.1 ms all the trapped particles are lost and only passing particles are confined. In the simulations without the TBM fields another 0.6% of the initial tritons is lost due to pitch-angle scattering resulting in a confined triton population fraction of 45.6% after 100 ms. When the TBM fields are present, significant losses occur after 0.1 ms. Those losses vanish at around 100 ms as can be seen in figure 13. There are two reasons for the vanishing loss rate: (i) the region of phase space where TBM-induced losses can come from is depleted and because the tritons are well above the critical energy, the depleted phase-space region at high triton energies is not refilled by pitch-angle scattering of low-energy tritons and (ii) the effect of slowing-down on the drift orbits is a shift of the orbits to the plasma center and hence co-going particles move away from the low-field side plasma edge where the TBM error fields are located. At the end of the simulation with the TBM fields only 30.6% of the original triton population is still in the plasma which is a reduction of 33% compared with the run without TBM fields present. This reduction compares very well with the observed value of 33 ± 5% as reported in [28].

5.3. RMP-induced beam-ion losses

RMPs can be beneficial for controlling ELM activity but they can have a detrimental effect on the fast-ion confinement, leading to increased losses and fast-ion redistribution. The SPIRAL code has been used in the planning phase of the

Figure 13. The confined triton fraction in the SPIRAL simulations as a function of time for a simulation without (green) and with the TBM fields present (red).

experiments to investigate at which RMP field amplitudes the effects of fast-ion loss and redistribution become observable in NSTX experiments.

For this study a standard NSTX high-density H-mode discharge was used. The fast-ion beam-birth profile was obtained from TRANSP/NUBEAM for beam source C with a tangency radius of 0.49 m (figure 14). Particle starting times were distributed uniformly between zero and 25 ms and all the 10 000 particles were then followed until they got lost or reached 25 ms. Slowing-down and pitch-angle scattering were included in the calculations where the energy slowing-down time varied from 23 ms in the core to less than 12 ms near the edge. When the vacuum-field amplitude of the non-resonant RMPs which were added perturbatively to the equilibrium magnetic fields was varied the beam ion losses increased as can be seen in figure 15. The RMP field is expressed as coil current where a coil current of 1.5 kA gives a radial magnetic perturbation of 3 mT at the LFS plasma edge as can be seen in figure 16(a). The RMP-induced losses are concentrated just below the mid-plane at the LFS. When we integrate the losses to the low-field side wall between the mid-plane and 0.5 m below the mid-plane as a function of toroidal angle we can clearly see that the loss pattern is modulated with the n = 3 RMP as is shown in figure 16(b) for the 1.5 kA calculation. Peak heat loads in that region of up to 200 kW m-2 for a 1MW beam are predicted when the RMPs are present while the heat loads are less than 20 kW m-2 for a 1 MW beam without the RMPs present. The minimum heat load that an infrared camera which is installed on NSTX can detect is less than 100 kW m-2 so it is concluded that RMP generated heat loads can be detected.

With the FIDA diagnostic the confined beam-ion population can be measured [15]. The interpretation of the FIDA signals is not straightforward because of its complicated weighting function in phase space. In the FIDA system on NSTX vertical viewing cords are used which favor low-pitch particles (\v\\/v\ < 0.6) in an energy range between 30 and 60 keV. For a quantitative interpretation of the FIDA signals, fast-ion distributions from amongst others TRANSP

NSTX pulse 142293

0.5 1.0 major radius [m]

Figure 14. NSTX plasma shape with the beam birth profile as obtained from TRANSP/NUBEAM for source C with a tangency radius of 0.49 m.

—i-1-1-1—

—i-1-1-1—

NSTX pulse 142293

—i-1-1-1—

CO CD .Q

_i_i_i_i_

_i_i_i_i_

_i_i_i_i_

0.0 0.5 1.0

RMP coil current [kA]

NSTX pulse 142293

90 180 toroidal angle |

1........1........1

with RMP 1

NSTX pulse 142293

no RMP

1........1 .

90 180 toroidal angle

Figure 15. RMP-induced beam ion losses as a function of the current through the RMP coils.

and SPIRAL can be used in theFIDASIM code [30] to calculate the FIDA response to those distributions. For the RMP loss study we have used the fast-ion distribution as calculated with SPIRAL (figure 17(a)) in the FIDASIM code to predict the

Figure 16. (a) The radial RMP field at the plasma boundary and (b) calculated beam-ion losses below the LFS mid-plane without RMP (blue) and with 1.5 kV RMP (red).

FIDA light that is observable as shown in figure 17(b). From this figure it can be seen that the RMPs decrease the confined beam-ion population sufficiently to detect changes in theFIDA light. Therefore, based on the SPIRAL simulations, it was concluded that enhanced RMP induced beam-ion losses can be detected experimentally on NSTX from both temperature measurements of the first wall for the lost ions and FIDA measurements for the confined fast ions.

6. Time-varying magnetic and electric fields

Simulating the interaction between electromagnetic waves and fast ions is of paramount importance to understand the fast-ion transport in tokamaks and therefore, a spectrum of electromagnetic waves can be included in SPIRAL calculations where we have made a distinction between ideal MHD waves that have no electrical field component in the direction of the local magnetic field and radio-frequency (RF) waves in the ICRF that are launched into the plasma for heating and/or current drive purposes. In general, RF waves have an electrical field component parallel to the magnetic field.

NSTX pulse 142293

major radius [m]

Figure 17. (a) FIDA signal without (blue) and with (red) RMPs as obtained with the FIDASIM code with the (b) the beam-ion distributions calculated with SPIRAL as input.

In SPIRAL the electromagnetic waves are taken as harmonic waves in the toroidal direction (see [31]):

EE(R, v, Z, t) = e(R, Z) exp(-i(ny + at)), B(R, y, Z, t) = b(R, Z) exp(-i(ny + at)) (20)

with e(R, Z) and b(R, Z) the complex fields that describe the poloidal distribution of the wave fields, n is the toroidal mode number and a is the mode frequency. The electric field EE and magnetic field B are related to each other by Maxwell's laws.

Ideal MHD codes such as NOVA [31,32] calculate amongst other quantities the plasma displacement, f, on a mesh inside the plasma. The plasma displacement can be used to obtain the electrical and magnetic wave fields for the mode from

EE = — x B = -iaf x B, dt f ,

-= -iaB = Vx E

B = Vx(f x B). (21)

Evaluating the EE from f and B can be done accurately. However, extracting B from equation (21) is numerically highly inaccurate because the curl operator requires the subtraction of large but almost equal terms. These terms have to be known with a very high precision to calculate B accurately. Ideal MHD codes usually give as an output the fluctuating magnetic field which is calculated accurately inside NOVA. We have used the plasma displacement together with magnetic fields as calculated by NOVA as input of the MHD fluctuations for SPIRAL.

The plasma displacement and fluctuating magnetic fields are given by NOVA on a discrete mesh and we have used again Chebyshev polynomial expansions to interpolate those quantities. The fluctuating electrical field is obtained from the plasma displacement and equilibrium magnetic field followed by a Chebyshev polynomial expansion:

E = -iaf x B = J2 J2eijTi(R)Tj(Z). (22) i j

Depending on the size of the structures in the MHD fields, 50100 basis functions are needed for each coordinate. The ideal MHD condition, EE ■ B = 0 is fulfilled automatically in this way and tests have shown that the fitted electrical fields satisfy this condition better than 10-4. The accuracy of the ideal MHD condition was improved further during the orbit calculations by adding tiny corrections to the fluctuating electrical field that forces the parallel electrical field to zero. This correction for E|| is usually less than 0.1 Vm-1 while E± is typically 103 V m-1 or higher for TAEs.

In a similar manner the radial, BR, and vertical, BZ, components of the magnetic field as calculated by NOVA are interpolated with Chebyshev polynomials and the toroidal field component, Bv, is calculated analytically from the divergence free condition imposed by Maxwell's equations and the Chebyshev coefficients of BBR and BBZ. The integration over the toroidal angle is analytically simple because of its harmonic dependence (equation (20)). Numerically, the procedure for obtaining the toroidal field component is robust because it does not depend on subtle cancellations of large terms as was the case with the curl operator. When following a particle through the plasma we have to evaluate now two complex fluctuating fields, EE(R, Z) and B(R, Z), which can be done efficiently at each (R, Z) location: the Chebyshev basis functions are evaluated once followed by the multiplication of the coefficients for each of the twelve field components. The complex fields are then multiplied by the toroidal and harmonic dependence, exp(-i(ny + at)) and the real part is used in the orbit-following calculations.

Fitting the fluctuating electrical and magnetic fields independently that were obtained in accordance with Maxwell's equations by NOVA does not guarantee that the fitted fields obey Maxwell's laws but where possible comparisons were made and the fitted fields were consistent. Moreover, It can be seen from equation (21) that the

1.0 1.5 2.0 1.0 1.5 2.0 1.0 1.5 2.0

major radius [m]

Figure 18. The radial (a), surface (b) and parallel (c) magnetic field components together with the radial (d) and surface (e) components of the electrical field of a RSAE. In (f) a single deuteron orbit. The full orbit is shown in black while the guiding center orbit is shown in gray. The radial contours are minor radius contours.

electrical field fluctuation amplitude is proportional to the mode frequency while the magnetic fluctuation amplitude is independent of the mode frequency. Hence, for low-frequency MHD activity the perturbed magnetic fields are dominant while at high frequencies the perturbed electrical fields become dominant. Therefore, small deviations of the independently fitted fields from the constraints of Maxwell's equations are tolerable because usually only the magnetic or the electrical field fluctuations are dominant. The alternative, evaluating the curl operator in equation (21), is less accurate as discussed earlier.

6.1. Single partice resonances

A reversed shear Alfven eigenmode as calculated with the NOVA code is shown in figure 18 where the radial (a), surface (b), and the parallel (c) components of the fluctuating magnetic field are shown together with the radial (d) and surface (e) components of the electrical field. The parallel electrical field is zero within the numerical accuracy of the calculation. In figure 18(f) a deuteron orbit is shown where the deuteron was followed for 1 ms. During this time interval the mode frequency increased linearly from 67.5 to 93.5 kHz. The

particle was launched on the mid-plane at a major radius of 1.95 m with an initial energy of 74 keV. After a number of poloidal passages the particle starts to resonate with the mode during a number of poloidal transits as can be seen from the velocity and energy traces in figure 19. It quickly gains energy (figure 19(c)) with a significant change in perpendicular velocity (figure 19(a)) while the parallel velocity is not affected (figure 19(b)) and its orbit shrinks by about 3 cm as can be seen in figure 19(d).

Resonant interactions between MHD modes and fast particles are important for fast-ion transport in plasmas where the linear resonances are given by

omhd = po + no^ (23)

with oMHD the frequency and n is the toroidal mode number of the mode, ot is the particle poloidal transit frequency, is the toroidal transit frequency and p is the bounce harmonic number [33-35]. At a given location, the transit frequencies depend on the energy and pitch. The resonances can be determined in energy-pitch phase space by scanning the energy and pitch at a fixed initial location as is shown in the red diamonds in figure 20 for the RSAE of figure 18 at 85 kHz where the poloidal and

0.0 0.2

0.4 0.6 time [ms]

0.8 1.0

0.0 0.2 0.4 0.6 0.8 1.0 time [ms]

Figure 19. Time traces of (a) the perpendicular and (b) parallel velocity, (c) the energy and (d) the major radius of a deuteron moving in the RSAE shown in figure 18.

^ 60 .y

0 -1.0

0 c o o

CD W O

0.0 pitch: v||/v

20 40 60 80 100 initial particle energy [keV]

Figure 20. Linear wave-particle resonances for the RSAE shown in figure 18 for particles launched at the mid-plane at a major radius of 1.95 m. The red diamonds were determined from the toroidal and poloidal transit frequencies while the blue curves are a measure for the energy exchange between the particle and the mode where an excursion of 0.1 in pitch corresponds to a change in energy of 1 keV.

toroidal transit frequencies were determined with the Fourier method as discussed above.

The resonance condition, however, does not give information on the strength of the interaction between the particle and the wave. This information can be obtained by studying the energy exchange between the mode and the particle. As a numerical experiment, we have launched particles at the mid-plane at a major radius of 1.95 m, scanned the pitch and energy, and recorded the minimum and maximum energy excursion that the particle makes during a 3 ms time interval due to the RSAE where we have used a very low RSAE amplitude with a maximum density fluctuation of n/n = 10-4

Figure 21. Energy exchange between a single deuteron and the RSAE shown in figure 18 as a function of the initial particle energy for a mode amplitude of (a) n/n = 1 x 10-4 and (b) n/n = 4 x 10-3. The particles were launched at the mid-plane at R = 1.95m and apitch of v^/v = 0.5.

to obtain only the linear resonances as given by equation (23). From figure 20 it can be seen that the locations where the mode exchanges energy with the RSAE coincide very well with the resonances as calculated from equation (23). A good coupling between the particle and the mode occurs for resonances with a pitch between -0.25 and 0.05 and 0.25 and 0.75. Outside those pitch ranges the interaction between wave and particle becomes weak.

When we increase the RSAE amplitude to an experimentally more realistic value of n/n = 4 ■ 10-3, the resonances broaden as can be seen in figure 21 and new fractional resonances appear due to the non-linear interaction between the particle drift orbit and the mode as explained

Figure 22. Time traces of the (a) measured TAE mode frequency and (b) measured amplitude of a burst of four TAEs in NSTX labeled with their toroidal mode number. The TAEs induce (c) fast-ion loss rate and (d) a redistribution of the fast ions at the mid-plane. The yellow line in (d) indicates the location of the magnetic axis. The (e) n = 2, (f ) n = 2, (g) n = 4 and (h) n = 6 TAE density fluctuation structures as calculated with NOVA at the maximum mode amplitude. Note the different color scales.

1.0 time [ms]

0.5 1.0 1.5

major radii

0.5 1.0 1.5 radius [m]

in [36]. The effect of the increased number of resonances and the broadening is that the resonances can overlap more easily, leading to chaotic particle orbits and enhanced fast-ion transport [36,37].

6.2. TAE-inducedfast-ion redistribution

TAEs appear regularly in NSTX as short (less than 2 ms as shown in figure 22(b) bursts when NBI heating is applied. Those modes reach large transient amplitudes and the mode frequency is changing rapidly as can be see in figure 22(a). The curves shown in figures 22(a) and (b) are polynomial fits to the measurements that are used in the SPIRAL simulations whereby the mode amplitude was determined from reflectometer measurements and the frequencies were measured with Mirnov coils. Those modes affect the fast-ion population and in the following we will show how the SPIRAL code can model the fast-ion redistribution in the presence of those TAE bursts.

In the time period preceding the TAE bursts the plasma was quiescent while 6 MW of NBI was injected from all three beam sources so that the fast-ion distribution at the onset of the TAE burst can be described well as slowing-down distribution. This slowing-down distribution was calculated with SPIRAL from the NBI-birth profile as obtained from the TRANSP/NUBEAM codes by loading 30 000 particles from all three energy components uniformly in time between 0 and 50 ms. Those particles were then followed with SPIRAL until 50 ms where slowing-down and pitch-angle scattering were

included. Particles that approach the thermal energy, in this case set at twice the local temperature, before they reached 50 ms were not followed any longer and marked as thermalized. The location of the particles in the parallel and perpendicular velocity phase space is shown in figure 23(a).

This distribution was then modified by the burst of TAEs after 2 ms as shown in figure 23(b) while particles were being added from the birth profile with the same rate as during the first 50 ms (i.e. the beams continued to fuel the plasma during the TAE bursts with the same rate as before). The mode amplitudes and frequencies varied continuously with time in the simulations according to curves shown in figures 22(a) and (b). It can be seen in figure 23(b) that some particles get accelerated in both the parallel and perpendicular direction to above the beam-injection energy and that a fraction of the particles, mainly with low pitch (vy/v) values get lost to the plasma boundary. The lost particle rate clearly coincides with the onset and growth of the modes as can be seen in figure 22(c) and in total 18% of the fast ions are expelled in 2 ms.

The modes have also a significant effect on the confined fast-ion distribution as is depicted in figure 24 where we show the evolution of the confined fast ions between -25 and 25 cm around the mid-plane. When the n = 2 TAEs amplitude increases after 0.2 ms and the n = 4 and 6 TAEs appear, the fast-ion profile collapses and broadens toward the high-field side. This is shown more clearly in figure 22(d) where the mid-plane distribution minus the distribution at the on-set of the TAEs is shown. The collapse of the peak is visible as a

J 2 £= CD CP

CD 1 CP

-10 12 3 parallel velocity [106 m/s]

TD £=

'WliiSm

parallel velocity [106 m/s]

Figure 23. The location of the particles in parallel/perpendicular phase space for (a) a beam-ion slowing-down distribution at the onset of the TAE burst and (b) after the TAE burst. The black dots represent confined particles while the orange dots represent particles lost to the last closed flux surface. The beam injection energy was 90 keV (edge of yellow circle) while the half and third energy components are at the boundaries of the purple and green circles. The red circle in the middle represents the thermal plasma and is drawn at 3 keV which is three times the central electron temperature.

—i-1-1-1—

—i-1-1-i—

—i-1-1-1—

—i-1-1-1—

passing

iad\usH

Figure 24. Evolution of the confined fast-ions between -25 and 25 cm around the mid-plane due to the TAEs as shown in figure 23.

trough (in blue) that forms at a major radius of 1.1 m while a ridge (in white) becomes visible around 0.8 m when the TAEs are excited. Such changes can be measured with the FIDA diagnostics.

The change in the distribution is a reflection of changes in the different particle classes as shown in figure 25. When the TAEs become active a fraction of particles on stagnation orbits which are located on the low-field side are changed into passing orbits which traverse the high-field side. This explains the broadening of the distribution. The particle population on stagnation orbits is also found to be sensitive to getting lost which contributes significantly to the flatting of the distribution as shown in figure 24. The fraction of trapped particles and particles on potato orbits is only slightly affected by the burst of TAE activity as can be seen in figure 25.

6.3. ICRFfields

In contrast to ideal MHD modes, externally launched RF fields for plasma heating and current drive have in general

stagnation trapped

1.0 time [ms]

Figure 25. Evolution of the confined fast ions between -25 and 25 cm around the mid-plane due to the TAEs as shown in figure 23.

an electrical field component parallel to the magnetic field direction. The full-wave code TORIC [38] can calculate the ICRF fields in tokamaks and its output routines were adapted to create an output file with RF fields that can be read directly by SPIRAL. In SPIRAL the TORIC fields are fitted in the same way as for the MHD modes but without the restriction of a vanishing parallel electrical field component. Because SPIRAL is a full-orbit code, the interaction between the gyro-motion of the particle and ICRF fields is included naturally. This is different from guiding center codes where a kick-operator has to be constructed to include the RF interaction [8,9]. In this section, we will show results from

_l_I_L

J_I_I_

—4.6

4.0 6.0 8.0 major radius [m]

ICRF simulation with SPIRAL for ITER

40 60 80 100 time [|os]

Figure 26. (a) The radial electric field for 20 MW of ICRF heating in ITER at 52.5 MHz where the fundamental 3He resonance is located near the plasma center (white). A trapped 500 keV deuterium orbit is indicated in purple together with its resonance at the high-field side (black). (b) The deuterium energy, (c) the gyro-frequency at the particle location and (d) the particle's parallel velocity for the depicted orbit.

two simulations in which RF fields are included. The first one is a study of ICRF absorption by NBI ions in ITER in which only the fundamental resonance layer exists in the plasma [39]. In the second study results are presented for HHFWs in NSTX where multiple overlapping resonances are present in the plasma [40].

In the TORIC calculations that were used below, the effects of the fast-ion currents generated on the dielectric tensor were calculated from a Maxwellian distribution. As it will be shown below, the fast-ion distribution is highly anisotropic. For an accurate comparison between the simulations and experiments the TORIC should include those non-Maxwellian distributions.

6.4. ICRF heating in ITER

One of the plasma heating schemes for ITER is a 20 MW ICRF system. This system can be tuned to 52.5 MHz which is the on-axis fundamental resonance frequency of 3He minority plasma species (figure 26(a)). The 3He should absorb most of the power but the 52.5 MHz deuterium resonance is also present in the plasma near a major radius of 5.2 m as shown in figure 26(a) and deuterium ions which form a slowing-down distribution from the 1 MeV beams can absorb some of the ICRF power thereby reducing the heating efficiency of the on-axis 3He minority heating. In figure 26(a) an orbit of a deeply trapped 500 keV deuterium particle is shown that passes through the resonance layer. The energy of this particle during one poloidal transit is shown in figure 26(b) together with the

499.8-

499.6-

time [ | s]

Figure 27. Details of the energy trace for the first crossing of the ICRF resonance layer as shown in figure 26(b).

local deuterium cyclotron frequency at the particle location (figure 26(c)) and its parallel velocity (figure 26(d)). From figure 26(b) it can be seen that the particle gets a kick in energy when it crosses the Doppler shifted resonance layer where the duration of one kick is about 1 \xs or 50 gyro-periods as can be seen in figure 27. Note also that away from the resonance a beat wave between the gyro-motion of the particle and the

Figure 28. (a) The radial electric field for 1 MW of HHFW heating in NSTX at 30 MHz. The different harmonic layers are shown in white and labeled with their harmonic number. The orbit in yellow is for a trapped 20 keV deuteron while its instantaneous guiding center is depicted in blue. (b) The deuterium energy, (c) the gyro-frequency at the particle location and (d) the particle's parallel velocity for the depicted orbit.

ICRF frequency can be seen. The Doppler effects on high-energy particles are significant; in this case the resonance is shifted from 52.5 to 50 MHz on the outer leg while resonance occurs at 54 MHz on the inner leg. From figure 26(a) it can be seen that the 50 and 54 MHz resonance layers are some 2030 cm away from the 52.5 MHz layer and therefore, it can be concluded that the ICRF resonance layers for energetic beam ions in ITER can be substantially broadened due to Doppler effects.

For a high 3He minority heating energy efficiency the deuterium beam ions should not absorb energy from the ICRF field. We have calculated the energy exchange between the ICRF waves and the beam ions by following 10 000 particles drawn from the beam-ion deuterium slowing distribution as obtained from TRANSP for 1 ms in SPIRAL with the 20 MW RF field as modeled with TORIC. From the energy gain of this ensemble the total absorbed power for the beam ions was calculated as 0.08 MW which is small compared with the 20 MW input power from ICRF and so it can be concluded that ICRF power absorption on beam ions in ITER is not important. The 0.08 MW compares well with the 0.06 MW as found with the TORIC code in which the fast-ion distribution was approximated by a Maxwellian distribution.

6.5. HHFWheating in NSTX

In contrast to conventional tokamaks, multiple resonances for HHFW heating are present in spherical tokamaks like NSTX as can be seen in figure 28(a) where all the resonances between

the second and the 11th harmonic are present. The drift orbits of the particles usually intersect more than one resonance and therefore, the particles can interact with multiple resonances. In figure 28(a) the orbit of a 20 keV deuteron is shown that was launched at the upper turning point on the 5th harmonic resonance layer. Despite crossing all of the harmonic layers between five and eleven, it only exchanges energy at the 5th harmonic layer (figures 28(b) and (c)) where it bounces and its parallel velocity is zero as can be seen in figure 28(d). The duration of the kick in this case is about 3 fs or 18 gyro-periods. At the other resonance layers the particle lingers for much shorter periods, often less than one or two gyro-periods, and therefore, it does not interact effectively with those layers.

Similar to the ITER case, the HHFW can couple to a population of fast ions that is created with NBI heating in NSTX and reduce the HHFW heating efficiency. We have used the SPIRAL code to estimate the amount of power that is absorbed by the beam ions. For this study we calculated the fast-ion slowing-down distribution from the beam deposition profile by loading 20 000 particles uniformly over a 50 ms time interval. Those particles were then followed to 50 ms without the RF fields present resulting in the slowing-down distribution shown in figure 29(a). This distribution was then followed from 50 to 51 ms in the presence of the 1 MW of HHFW power. During this interval beam ions from the birth profile were also loaded with the same rate as in the first 50 ms to simulate the continued beam fueling. The resulting distribution in parallel/perpendicular phase space is shown in figure 29(b) where the black dots represent particles that were confined

parallel velocity [106 m/s] parallel velocity [106 m/s]

Figure 29. The location of the particles in parallel/perpendicular phase space for (a) a beam-ion slowing down distribution just at the start of HHFW injection and (b) after applying 1 MW of HHFW heating for 1 ms. The black dots represent confined particles while the orange dots represent particles lost to the first wall. The beam injection energy was 90 keV (edge of yellow circle) while the half and third energy components are at the boundaries of the purple and green circles. The red circle in the middle represents the thermal plasma and is drawn at 3 keV which is twice the central temperature.

perpendicular energy [keV] parallel energy [keV]

Figure 30. (a) The perpendicular and (b) parallel energy spectra of all the particles before (red) and after (green) 1 ms HHFW heating. In blue the spectrum of the confined particles at 1 ms is shown.

at 51 ms while the particles that were lost during the 1 ms HHFW heating phase are marked with orange dots. From this figure it can be seen that the particles only gain perpendicular velocity from the HHFW in contrast to interactions with MHD activity where both parallel and perpendicular velocity changed (figure 23).

The parallel and perpendicular energy distribution spectra (figure 30) show the formation of a high-energy tail in the perpendicular energy (figure 30(a)). The majority of the accelerated ions are lost during the 1 ms HHFW heating interval while the parallel confined energy distribution (figure 30(b)) is hardly changed by the HHFW. The lost ions, however, have gained parallel energy which can be understood from their orbits. The parallel energy of trapped particles varies between zero at the bounce points to a maximum on its outer leg as can be derived from the velocity trace shown in figure 28(d). The HHFW increases the perpendicular particle energy near

the bounce point. This increases the size of the drift orbit and the maximum parallel energy of the particle at its largest major radius. Those accelerated particles then get lost at the LFS near the mid-plane with an increase in parallel energy that was obtained from HHFW at the bounce point.

From the 1 MW of HHFW power that was used in this simulation 21 kW (or 2.1% of the injected HHFW power) was absorbed by the fast ions and 16 kW of that power was lost to the wall. Although a small amount of HHFW power is carried away by the accelerated fast ions, the impact of the waves is very severe for the fast-ion population. In the 1 ms when the HHFW are present 13% of the fast particles are lost. In a separate run where NBI injection into a HHFW heated plasma was simulated for 5 ms it was found that up to 53% of the injected particles were lost while in the absence of the HHFW power only 13% of the ions were lost as first-orbit losses. Therefore, HHFW heating in combination with NBI

700 600

10 15 20 25 Larmor radius [cm]

Figure 31. Distribution of the Larmor radii for the confined particles at the start (red) and at the end of HHFW heating (blue). The Larmor radius distribution of the confined and lost particles are indicated in green.

150 100 50 0

1 1 1 1 1 1 1 1 1 1 i i i i i i i i

■Ifff> t V V- Mi ' TÜí"

- , ......... i

0 20 40 60 80 100 energy at HHFW start [keV]

Figure 32. Energy gain as function of the particle energy at the start of the 1 ms HHFW heating for an NBI slowing-down distribution as shown in figure 29. The black marks represent particles that are confined after 1 ms while the orange dots represent the lost particles.

heating seems to lead to a decreased plasma performance in NSTX.

The perpendicular energy gain is also reflected in the increase in Larmor radii as shown in figure 31. The separation between the different HHFW harmonics (figure 28(a)), about 16 cm in this simulation, is comparable or smaller that the gyro-radii of the high-energy particles and therefore, particles can interact with two resonance layers at one gyro-orbit. The effects of this harmonic overlap is currently being investigated theoretically [41].

In order to model the interaction between HHFW and NBI accurately, the effects of the fast-ion currents generated by the NBI and HHFW on the dielectric tensor must be included in the calculation of the HHFW fields. In the TORIC code the dielectric tensor is currently calculated from a Maxwellian distribution. For an accurate comparison between the simulations and experiments the TORIC code is being modified to include non-Maxwellian distributions in the calculation of the dielectric tensor [42]. Because the RF waves modify the fast ion distribution which in its turn modifies the dielectric tensor and hence the TF fields in the plasma, an iterative procedure between SPIRAL and TORIC is being developed to study the HHFW-induced modifications of the fast-ion distribution and its feedback to the RF fields in the plasma for detailed comparison between simulations and experiments.

The HHFW couples most efficiently to the high-energy particles from the slowing-down distribution and not to the bulk thermal plasma as can be seen in figure 32 where particles with initial energies above 10keV get accelerated significantly. The thermal plasma temperature in the center used in the simulations was 1.5 keV while the distribution above 10 keV is part of the non-Maxwellian beam-ion slowing-down distribution. This again highlights the point that non-Maxwellian distribution functions should be taken into account in the RF-field calculations for a detailed comparison with experiments.

7. Conclusion

In this paper, we have given a description of the numerical methods used in the full-particle-orbit following SPIRAL code together with a number of physics studies performed with the code to illustrate its capabilities. Full-orbit effects are important for fast-ion experiments in low-field machines such as NSTX. The SPIRAL code is written as a tool to be used in the analysis and interpretation of fast particle experiments and it is also useful in the planning phase of experiments to investigate their feasibility.

Apart from the toroidally symmetric equilibrium magnetic fields, time independent toroidal ripple fields, electric fields arising from the plasma rotation, and time dependent magnetic electric fields from MHD modes and ICRF antennas can be included in the SPIRAL code. Pitch-angle scattering in 3D and slowing-down effects, which are important for fast-ion transport studies are also included in the code. Realistic first walls for different machines can be selected in the code so that particle losses and heat loads in the first wall can be studied accurately. A large number of variables is written in the output file (appendix B) for access with a post-processor code so that the results can be studied in detail and routines are available to explore the output further with the Interactive Data Language (IDL) [43].

In this paper, we have shown some results of studies performed with the SPIRAL code where full-orbit effects play a role. The exact location of hot spots created by fast-ion losses depends on the gyro-radius, especially in combination with strong magnetic field gradients which are present close to the wall in the TBM experiments in DIII-D. Finite gyro-orbit effects are also dominant for the calculation of the 1 MeV tritium burn-up experiments in the presence of the TBM fields.

Fast-ion gyro-radii in NSTX are large because of its low magnetic field while the typical scale of TAEs are on the order or smaller than those gyro-radii which justified a full-orbit

treatment for the interaction between TAEs and fast ions. It was found that a short burst of TAEs can redistribute the fast ions and induce significant losses.

The interaction between ICRF and/or HHFW and fast ions depends solely on the gyro-motion of the fast ions. We have shown for an ITER plasma that the wave-particle interaction occurs naturally in the SPIRAL calculations while in a preliminary study for NSTX we found that the HHFWs couple efficiently to the energetic beam ions. The beam ions got accelerated where their gyro- and drift-orbits increase so that they do not fit into the machine any more and get lost to the wall.

Although, where possible, comparisons between full-orbit and guiding-center calculations were made, an exhaustive study between the two calculations is beyond the scope of this paper. For unambiguously studying the differences between the two approaches one has to make sure that the differences do not arise from small differences in the magnetic and/or electrical fields in the different codes.

In a similar way, the derivative of a Chebyshev polynomial can be written as a sum of lower order polynomials:

— Ti(x) = 2i y^ Tj(x) - iT0(x) i is odd, dx ^

j =0 step2

is even.

Ti(x) = 2i ^^ Tj(x) i is

j = 1 step2

When an arbitrary function is written as a finite sum of Chebyshev polynomials:

f(x) = ^2 a Ti(x)

the integral and derivative can be obtained analytically as finite Chebyshev polynomial by rearranging the coefficients, ai, with the proper weights from the equations (A1)or( A2). When the function depends on two independent variables, f (x, y), the expansion coefficients form a matrix:

Acknowledgments

This work was supported by the US Department of Energy under DE-AC02-09CH11466, SC-G903402 and DE-FC02-04-ER54698.

f(x,y) = aijTi(x)Tj(y)

i=0 j=0

and integrals and derivatives can be obtained in a similar way as for the one variable case.

Appendix A. Chebyshev polynomials

Appendix B. Code output

The Chebyshev polynomials are defined on the interval [-1.0, 1.0] by the following recursion relation:

Ti+i(x) = 2xTi(x) - Ti-i(x)

with T0 = 1.0 and T1 = x and they form an orthogonal set on their interval. Their range is bounded by [-1.0, 1.0] and together with the fact that the extrema are distributed uniformly the Chebyshev polynomials are well suited for numerical work. When a function is expressed as a sum of Chebyshev polynomials:

f(x) = OiTi(x),

where a simple coordinate transformation is made to map the finite range of the function onto the interval [-1.0, 1.0] where the Chebyshev polynomials are valid. When the function is needed up to a certain numerical accuracy, e, the sum can be truncated from where the coefficients, ai, are less than e. The uniform spread of the extrema then guaranties that the function is accurate up to e on its whole range.

The integral of a Chebyshev polynomial can be expressed as a the sum of two Chebyshev polynomials:

T¿(x) dx =

-Ti+1(x)--,

2(i + 1) 2(i - 1)

4 T2(x)

i = 1,

|(x) i = 1,

where one of the polynomials is one order higher and the other one order lower than the original polynomial order.

When individual orbits are calculated with SPIRAL as shown in section 2 the position and velocity of the particle are written in the output file as a function of time together with the magnetic electric and fields and the poloidal flux at the particle location which is sufficient to obtain various other quantities that are important for the fast particle dynamics. When ensembles of fast ions are calculated it is impossible to store information for each orbit individually because of disk space constraints. In this case only the information at the start and end times is stored which consists of the particle location and velocity, the magnetic field and poloidal flux at the particle location and at the guiding center as calculated with equation (4) and the change in particle energy due to slowing down and interactions with the electrical fields. At the start and end time the poloidal and toroidal transit frequencies are calculated and the particle orbits are classified which is also written on the output file. This information is sufficient to calculate a large number of other quantities such as magnetic moments, Larmor radii, toroidal angular momenta, particle pitches, etc. which can be accessed with the post-processor interface for IDL.

References

[1] Northrop T G 1961 Ann. Phys. 15 79

[2] Littlejohn R G 1981 Phys. Fluids 24 1730

[3] Littlejohn R G 1983 J. Plasma Phys. 29 111

[4] White R B 2006 The theory ofToroidally Confined Plasmas

revised 2nd edn (London: Imperial College Press) ISBN 1-86094-639-9

[5] Gardner C S 1966 Phys. Fluids 9 1997

[6] Weyssow B and Balescu R 1986 J. Plasma Phys. 35 449

[7] Hauff T, Pueschel M J, Dannert T and Jenko F 2009 Phys. Rev.

Lett. 102 075004

[8] Choi M, Chan V S, Pinsker R I, Chui S C, Heidbrink W W

2005 Phys. Plasmas 12 072505

[9] Choi M, Chan V S, Berry L A, Jaeger E F, Green D, Bonoli P,

Wright J and the SciDAC Team 2009 Phys. Plasmas 16 052513

[10] Kernighan B and Ritchie D M 1988 The C Programming

Language 2nd edn (Englewood Cliffs, NJ: Prentice Hall) ISBN 0-13-110362-8

[11] The Numerical Algorithms Group 2012 http://www.nag.co.uk/

[12] Gropp W, Lusk E and Skejellum A 1999 Using MPI: Portable

Parallel Programming with the Message-Passing Interface 2nd edn (Cambridge, MA: The MIT Press) ISBN 0-262-57134-X

[13] Boozer A H and Kuo-Petravic G 1981 Phys. Fluids

24 851

[14] Perkins R J et al 2012 Phys. Rev. Lett. 109 045001

[15] Heidbrink W W, Burrell K H, Luo Y, Pablant N A

and Ruskov E 2004 Plasma Phys. Control. Fusion 46 1855

[16] Lao L L, Ferron J R, Groebner R J, Howl W, St John H,

Strait E J and Taylor T S 1990 Nucl. Fusion 30 1035

[17] Mason J C and Handscomb D 2003 Chebyshev Polynomials

(Boca Raton, FL: Chapman and Hall/CRC) ISBN 0-8493-0355-9

[18] Fu G Y, Park W, Strauss H R, Breslau J, Chen J, Jardin S and

Sugiyama L E 2006 Phys. Plasmas 13 052517

[19] Cordey J G and Core W G F 1974 Phys. Fluids 17 1626

[20] Chen F F 1983 Introduction to Plasma Physics revised 2nd

edn (New York: Plenum Press) ISBN 0-306-41332-9

[21] Huba J D 2009 NRL plasma formulary, Naval Research

Laboratory Washington, DC

[22] Budny R V et al 1995 Nucl. Fusion 35 1497

[23] Pankin A, McCune D, Andre R, Bateman G and Kritz A 2004

Comput. Phys. Commun. 159 157

[24] Bortolon A et al 2013 Redistribution of Fast Ions by Low

Frequency MHD Instabilities in NSTX H-mode Plasmas, in preparation

[25] Goldston R J, White R B and Boozer A H 1981 Phys. Rev.

Lett. 47 647

[26] Tobita K, Nakayama T, Konovalov S V and Sato M 2003

Plasma Phys. Control. Fusion 45 133

[27] Giancarli L et al 2010 Fusion Eng. Des. 85 1829

[28] Schaffer M J et al 2011 Nucl. Fusion 51 103028

[29] Kramer G J et al 2011 Nucl. Fusion 51 103029

[30] Heidbrink W W, Liu D, Luo Y, Ruskov E and Geiger B 2011

Commun. Comput. Phys. 10 716

[31] Cheng C Z 1992 Phys. Rep. 211 1

[32] Cheng C Z and Change M 1987 J. Comp. Phys. 71 124

[33] Pinches S D, Kiptily V G, Sharapov S E, Darrow D S,

Eriksson L-G, Fahrbach H-U, Garcia-Murioz M, Reich M, Strumberger E, Werner A and the ASDEX Upgrade Team and JET-EFDA Contributors 2006 Nucl. Fusion 46 S904

[34] Heidbrink W W 2008 Phys. Plasmas 15 055501s

[35] García-Muñoz M et al and the ASDEX Upgrade Team 2009

Nucl. Fusion 49 085014

[36] Kramer G J, Chen L, Fisher R K, Heidbrink W W, Nazikian R,

Pace D C and Van Zeeland M A 2012 Phys. Rev. Lett. 109 035003

[37] García-Munoz M et al and ASDEX Upgrade Team 2010 Phys.

Rev. Lett. 104 185002

[38] Brambilla M 2002 Plasma Phys. Control. Fusion 44 2423

[39] Budny R V et al and Members of the ITBP-IOS 2012 Nucl.

Fusion 52 023023

[40] Phillips CK et al and the NSTX Team 2009 Nucl. Fusion

49 075015

[41] Burby J W, Kramer G J, Phillips C K, Valeo E J 2011 AIP

Conf. Proc. 1406 345

[42] Bertelli N 2012 private communication

[43] Interactive Data Language home page:

http://www.exelisvis.com/language/en-us/ productsservices/idl.aspx