Simulations of beam ion transport during tearing modes in the DIII-D tokamak

Large coherent MHD modes are observed to reduce the neutral beam current drive efficiency and 2.5 ,MeV neutron emission in DIII-D by as much as ~65%. These modes result in large (width w≲20 cm for minor radius a≈60 cm), stationary, single helicity magnetic islands, which might cause anomalous deuterium beam ion losses through orbit stochasticity. An analytic estimate predicts that co-going, passing deuterons with E≳40 keV become stochastic at island widths comparable to those in the experiment. A Hamiltonian guiding centre code is used to follow energetic particle trajectories with the tearing mode modelled as a radially extended, single helicity perturbation. In the simulations, the lost neutral beam current drive and neutron emission are 35% and 40%, respectively, which is consistent with the measured reductions of 40±14% and 40±10%. Several features of the lost particle distribution indicate that orbit stochasticity is the loss mechanism in the simulations and strongly suggest that the same mechanism is responsible for the losses observed in the experiment.


Introduction
Understanding fast ion confinement in the presence of perturbing helical fields is important for predicting the behaviour of energetic ions in future devices. Perturbations due to coherent magnetohydrodynamic (MHD) modes can lead to unwanted or unpredictable fast ion transport. Such anomalous losses would degrade the plasma performance and may damage the vessel wall and plasma facing components. In a reactor, early losses of energetic ions could prevent ignition while alpha particle losses would terminate the discharge during the burning phase.
Coherent MHD instabilities can affect fast ion confinement in several ways [1]. The strongest interaction is with fast ion driven instabilities. For these modes, the frequency of the observed mode is associated with some characteristic frequency of the fast ion motion and the mode gains energy at the expense of the fast ion population. Because the propagating wave and the fast ion orbit resonate, some ions stay in phase with the wave's electric field, resulting in large spatial transport.
Slowly rotating modes that are driven unstable by the background plasma are conceptually different. In this case, temporal resonances between the orbit and the mode are unimportant. Nevertheless, the static helical magnetic field perturbations associated with globally extended MHD modes can still degrade fast ion confinement [2][3][4][5]. The n = 0 orbit shift (where n is the toroidal mode number) of the fast particle motion couples to the added motion from the helical perturbations, resulting in drift islands in the particle's phase space that can overlap and cause the particle's motion to become stochastic. In this article, this transport mechanism is referred to as 'orbit stochasticity'. Orbit stochasticity is quite different from parallel transport along ergodic field lines. For ergodic field lines, finite gyroradius effects reduce the transport by effectively averaging over the fluctuations [6], while for orbit stochasticity, the finite gyroradius is responsible for anomalous transport. In theory, then, orbit stochasticity can cause or enhance anomalous losses of fast ions without affecting the confinement of the bulk ions. More importantly, fast particles on passing drift orbits, which are normally thought of as well confined, can be affected by orbit stochasticity.
Even though this transport mechanism is theoretically feasible and is suspected to have caused or enhanced fast ion losses in tokamaks during a variety of experiments [1], detailed experimental verification that orbit stochasticity causes fast ion losses is sparse. Zweben et al suggested that orbit stochasticity contributed to alpha particle losses on TFTR [7]. In that study, alpha particle loss measurements taken during a DT discharge with coherent MHD activity were compared with results from numerical simulations using a guiding centre code. Simulations for a range of mode strengths and MHD activity were consistent with orbit stochasticity, but problems in categorizing the MHD activity, uncertainties in modelling the vacuum field (which propagated into the loss measurements) and the difficulty of modelling the perturbation field made the agreement qualitative. A similar report appears in ref. [8]. In a study using externally imposed helical fields on DIII-D, no impact on fusion product confinement was observed, consistent with theoretical estimates [9]. In the Auburn torsatron, the effects of rotating, externally imposed magnetic perturbations on ion confinement were consistent with the predictions of intrinsic orbit stochasticity theory [10].
The comparison between theory and experiment described in this article is the most detailed experimental verification of orbit stochasticity as a mechanism for beam ion loss in a tokamak performed to date. During experiments at DIII-D, large, stationary MHD modes caused reductions in both the 2.5 MeV neutron production ( 65%) and the central noninductive current ( 55%) [11]. These data indicated losses of deuterium beam ions from the core without the confinement of the bulk species being affected. An analytic estimate shows that orbit stochasticity was important for beam ions in this discharge. Numerical simulations using a realistic equilibrium and initial beam ion distribution show quantitative agreement between neutron and current drive losses for magnetic island sizes comparable to those seen in the experiment. Features of the lost particle distribution function are consistent with orbit stochasticity as the loss mechanism in the simulations, implying that this same mechanism is responsible for fast ion losses seen in the experiment.
This article is organized as follows. Details of the experiment, including the loss measurements, the measurements of the magnetic islands and a discussion of experimental errors, are given in section 2. In section 3, an analytic estimate of the island overlap criterion is used to test if intrinsic orbit stochasticity is important for deuterium beam ions in this discharge. The details and results from the numerical simulations are given in section 4. Finally the implications and conclusions from this work are given in section 5.

Experimental data
Discharges from the neutral beam current drive (NBCD) experiment [11] with strong MHD activity consistently had much lower central non-inductive current drive and 2.5 MeV neutron emission than expected theoretically. For this study, a single discharge (shot 89384) from that experiment is considered. The discharge is a double null diverted, H mode deuterium plasma with 1.7 T toroidal field, a plasma current of I p 0.55 MA, a major radius of R 0 1.86 m at the magnetic axis, an elongation of 2.0 and a triangularity of 0.8. Near tangential (tangency radii R tan = 0.76 m and 1.15 m) deuterium beams inject ∼75 keV neutrals in the direction of the plasma current. The electron density is measured by Thomson scattering [12] and CO 2 interferometer [13] diagnostics, the electron temperature by Thomson scattering and electron cyclotron emission (ECE) [14], and the carbon ion temperature and density by the charge exchange recombination (CER) diagnostic [15]. (Carbon is the dominant impurity in this discharge.) The deuterium density and temperature are inferred from the electron and carbon data using the TRANSP code [16]. Representative profiles are shown in figure 1.
The goal in shot 89384 was to ramp up the plasma current using NBCD and bootstrap current. Initially, a 500 kA ohmic target plasma is created. At 1.2 s, the loop voltage V , or toroidal EMF, is fixed to limit the inductive current as the beam power is slowly increased, to build up the non-inductive current. Shown in figure 2 are the time evolution of the beam power, Mirnov activity and 2.5 MeV neutron emission for the discharge of interest, as well as data from a similar discharge with lower beam power, shot 89389. As the beam power and β (the ratio of the volume averaged plasma pressure to the magnetic field pressure) are increased, a large MHD mode with toroidal mode number n = 1 is detected by the Mirnov array at ∼1.95 s with a frequency of ∼21 kHz. No clear poloidal mode     The diagnostic has a sightline along the midplane and measures emissions along this chord. Channel 1 is located at the outboard midplane, and the channels are separated by 2 cm along the midplane chord. Phase shifts of 180˚and high radiated power are due to a large magnetic island (shaded region) centred at the q = 2 surface with width w = 17 ± 2 cm. The radial location of the q = 2 surface and the magnetic axis from a numerical equilibrium are noted in the figure. The open circles represent measurements whose coherency is too low to be statistically significant.
numbers are evident from analysis of the poloidal Mirnov array data. No high frequency MHD activity (f 50 kHz in the lab frame), associated with beam driven instabilities, is evident during this phase of the discharge. Between 1.5 and 2.5 s, the energy confinement time τ E stays level at ∼52 ms but starts to drop at 2.5 s and continues down to ∼44 ms at around 3.5 s.
Phase shifts of 180˚and relatively high power in the 19.8 kHz fluctuations from the ECE at 2.3 s (figure 3) are evidence of a large island chain about the q = 2 surface. (The dominant n = 1 mode from Mirnov data has a frequency of 19.8 kHz at this time.) This slice was selected since it has the clearest identification of the island structure, but the phase shifts near the q = 2 surface are evident continuously up until 2.8 s. (After 2.8 s, the character of the MHD activity changes.) At 2.3 s, the island width is w = 17 ± 2 cm along the outer midplane, compared with a nominal minor radius a = 59.9 cm. CER [15] indicates that the plasma is rotating at ∼21 kHz at this time, so these islands are stationary in the plasma frame. The stationary island structures indicated by the ECE at this frequency probably represent the same mode detected by the Mirnov array, so the poloidal mode number inferred from the ECE data is m = 2.
An independent measurement of the island width can be extracted from the magnetics data. The Mirnov signal from a probe at the vessel has a peak to peak amplitude of δB θ,wall 3.6 × 10 −4 T during the burst at 2.3 s. Using a multipole expansion about the magnetic axis, the magnitude of the radial component of the perturbation can be approximated by δB r 1 2 |B θ | wall (r probe /r peak ) m+1 , where r probe is the distance between the magnetic axis and the magnetic probe and r peak ρ peak a is the distance between the magnetic axis and the minor radius at which the mth poloidal harmonic peaks. (Note that ρ peak is the radial location of the peak of the mth poloidal harmonic in terms of the square root of the toroidal flux, ρ ∼ √ ψ.) The peak value of the radial component of the (2,1) mode is δB r /B = (4.1 ± 1.6) × 10 −3 , with the uncertainty dominated by the assumption that the mode's radial structure can be described by a multipole expansion. The saturated island width can be approximated using this value and equation (5) in ref. [17], resulting in a width w = 20 ± 10 cm. The uncertainty associated with this measurement is large, so at best it serves to check the consistency of the ECE data.
Low n, slowly growing MHD activity resulting in stationary magnetic islands at the q = m/n surface is the signature of a tearing mode. However, using a calculation of [18], the plasma is found to be conventionally stable for low (m, n) MHD modes, even considering transients in β [11]. ( is the usual resistive MHD stability parameter [19].) Most likely these tearing modes are destabilized by the pressure gradient modification of the neoclassical bootstrap current [20]. Assuming the existence of seed islands at some rational surface, the plasma pressure gradient reinforces the bootstrap current at the island location and acts as a free energy source for the tearing mode. The relatively high electron temperature (∼3 keV) and low density (∼5×10 13 cm −3 ) make neoclassical destabilization even more effective. (Similar pressure gradient driven tearing modes were observed previously on DIII-D [21] and TFTR [17].) The NBCD and neutron emission from this discharge are less than theoretically expected. The non-inductive current profile is given in figure 4, along with predictions from the TRANSP [16] and ONETWO [22] codes. (The inductive and non-inductive parts of the current density are separated by a technique which calculates the parallel electric field from time dependent equilibrium reconstructions and assuming neoclassical resistivity [23].) The bootstrap current is a significant contributor to the non-inductive current but the NBCD dominates in the core. The difference between theory and experiment for the NBCD current drive inside ρ 0.30, using a calculation for the bootstrap current, is 40 ± 14% for the TRANSP prediction and 43 ± 14% for the ONETWO prediction. The quoted uncertainty is due solely to the random error in the measured non-inductive current profile; not included are systematic uncertainties in the theoretical predictions. A significant reduction in the expected 2.5 MeV neutron signal starts when the mode turns on at 1.95 s. During the burst at 2.3 s, the neutron signal is ∼31% below a prediction from a 0-D code [24] that uses central plasma parameters to calculate the classical beam-plasma and thermonuclear neutron rate, as shown in figure 2. The signal is ∼40% below the prediction from the TRANSP code, which uses full experimental profiles and a numerical equilibrium from EFIT [25] to calculate a total neutron emission rate that also includes beam-beam reactions. TRANSP also finds that beam-beam reactions are significant for this discharge (∼30% of the total, with the thermonuclear contribution negligible), which explains the larger discrepancy of the data with the TRANSP prediction than the 0-D code prediction. For MHD quiescent discharges during the NBCD experiment, both the neutron emission and the non-inductive current profile compare well with the theoretical predictions. The difference between the theoretical predictions and the measured neutron rate is indicative of beam ion losses during the MHD activity. Using an ensemble of discharges with similar amplitude tearing modes to estimate the uncertainty, the reduction in neutron rate relative to the classical prediction is 40 ± 10%.
The plasma equilibrium is generated from a fit to the magnetics, motional Stark effect (MSE) [26] and kinetic data using the equilibrium fitting code EFIT. The MSE data used to constrain the equilibrium did not correct for the radial electric field E r , since those data were unavailable for this discharge. (E r has been shown to affect the MSE measurements for high performance discharges [27].) However, for these plasmas the E r correction is expected to be negligible, causing only a shift on the order of a centimetre in the location of the q = 2 surface. A slightly different q profile would not alter the efficacy of orbit stochasticity in the simulations (section 4.2). As discussed below, uncertainties in other quantities are expected to be much larger.
The uncertainties in the fast ion loss measurements and to some extent the measurement of the island width are the dominant sources of error in the comparison between experiment and theory. The systematic errors associated with the simulations discussed in section 4 are not expected to be nearly as large.
Consider first the uncertainty in the fast ion loss measurements. No direct measurements of the fast ion losses are available for this discharge. The absolute amount of fast ion losses is determined indirectly by obtaining a prediction from a calculation based on experimental data (in this case TRANSP) and then comparing this prediction with the experimental measurements (the current drive and neutron losses). Consequently there are uncertainties associated with both the calculation of the expected value and the measurements themselves which propagate into the overall current drive and neutron losses. In the case of the total current drive losses, the error bars on the measured noninductive current profile alone result in a ±14% uncertainty in the total expected current; the systematic errors in calculating the expected profiles are in addition to this amount. The neutron losses fare better. By taking an ensemble of discharges from the NBCD campaign with anomalous neutron losses and considering the scatter in the ratio of the measured to the predicted neutron emission, the measured neutron loss is 40 ± 10%, with the uncertainty including both systematic errors in the modelling and uncertainties in the measurements. The uncertainty in the island width is fortunately much smaller. The tearing mode results in such large magnetic islands that the signal to noise ratio in the ECE fluctuation data is excellent, as the islands span about eight or nine channels along the midplane. (These values translate to δρ = 0.26 ± 0.015 in flux co-ordinates.) The edge magnetics measurement of the perturbation, even with its large error bars, serves as an independent check that the width from the ECE measurements is reasonable.
The current drive reduction, neutron reduction and island width measurements are summarized in table 1, along with the associated uncertainties.

Estimate of island overlap
Using some experimental data, the importance of intrinsic orbit stochasticity in this shot can be estimated. An analytic approximation of the island overlap criterion for co-going, passing particles (the pitch λ ≡ v /v = 1 for such particles) in the presence of a single helicity island can be obtained from the q profile (figure 5) and the ECE measurements ( figure 3). The (2,1) island chain is centred on the q = 2 flux surface, corresponding to a radial position of ρ q=2 0.22, with a width δρ 0.26. (Technically, the island chain is offset towards the outer midplane from the q = 2 surface and is centred on the q = 2 drift surface, but in this estimate the offsets for q = 2 and q = 3 drift surfaces are assumed comparable, see equation (1).) The theory predicts that the perturbation results in a primary island in the particle's phase space at the q = 2 surface and sideband islands at the q = 1 and q = 3 surfaces  (ρ q=3 0.39). Since q > 1 for this discharge, only the q = 3 islands are considered. The widths of these drift islands in the particle's phase space are related to the magnetic island widths through δρ l δρ √ |G l (E, ε)|, where l denotes the island harmonic label (l = 0 is the primary island, l = ±1, ±2, . . . , are sidebands) while G l is the 'coupling coefficient' and is an implicit function of energy E and inverse aspect ratio ε [5]. The coupling coefficients are normalized such that |G l | 1 and the phase space islands will never exceed the width of the magnetic island. Using coupling coefficients calculated in ref. [5], G 0 0.95 and G 1 0.25 for the case of ∼40 keV deuterons in a large aspect ratio tokamak.
For the primary and sideband islands to overlap, the drift island half-widths must be greater than the distance between the centres of the drift islands, so the Chirikov island overlap criterion can be stated as The sum of the drift island half-widths comes out to ∼0.19 and the distance between flux surfaces is ∼0.20, meaning the overlap criterion is satisfied for ∼40 keV passing deuterons in this discharge. Higher energy particles have a lower stochastic threshold and most of the beam ions are on co-going, passing orbits, so this calculation indicates that orbit stochasticity is a good candidate for the loss mechanism in the experiment.

Methods and approximations
To study this transport mechanism numerically, the Hamiltonian guiding centre code ORBIT [28] is used to follow beam ion trajectories in a numerical equilibrium with static magnetic islands to model the effect of the tearing mode. The goal in the simulations is to obtain a 'steady state' distribution of fast particles from which current density and neutron emissivity profiles can be deduced. Initially, a 'birth' population with a single energy is created, very similar to a short beam pulse or 'blip'. ORBIT is used to follow fast ion trajectories, slowing down and pitch angle scattering in a numerical equilibrium. The position, pitch and energy of the confined population of ions are recorded at periodic time intervals, resulting in a series of 'snapshots' of the beam ion population. Then a steady state distribution of ions is assembled by combining these snapshots of the birth distribution as it evolved in time. To model the tearing mode, a zero frequency, radially extended, single helicity perturbation can be introduced in the simulations and allowed to act on the evolving birth distribution over time. The Monte Carlo birth particle population used in the simulation (figure 6) is obtained from TRANSP calculations of the beam deposition profile. The particles are assumed born along the outer midplane, θ = 0. Note that, in the TRANSP run, the particle population is calculated in the lab frame while the simulations are performed in the plasma frame. However, since the plasma rotation is <10% of the transit frequency of the beam ions, the implied Doppler shift in the particle velocities is negligible.
Both the current drive and the neutron measurements are sensitive to the energetic portion of the beam ion distribution function, so the simulations are run for a total of 18.2 ms, or 4200 transits for a centrally born 75 keV passing beam ion. (The transit time for such a particle is τ trans = const 4.3 µs.) After 18.2 ms, the probability that a typical central beam ion will produce a fusion reaction has dropped to <1/e of its initial value [29], and the probability rapidly decreases as the ion further decelerates.
Pitch angle scattering and energy slowing down are the only collisional effects included in the simulations. Pitch angle scattering is modelled using an energy conserving Monte Carlo operator [30], while energy slowing down is represented simply as an incremental change in energy over a time step. The pitch angle scattering collision frequency ν d and the energy loss rate ν E (figure 7) are calculated for every particle using interpolated experimental data (figure 1) in the standard Coulomb collision formulas [31]. The equations of motion are solved at fixed intervals of t = τ trans /30; the collision operators are applied directly afterward. Particles are checked to ensure that any energy change is under 20%, i.e. E/E < 0.2.
The position, pitch and energy of the confined ions are recorded at intervals of 100τ trans . Particles are considered lost when they cross the plasma edge in the numerical equilibrium. (Because drift orbits that cross the separatrix do not necessarily intersect a wall, this boundary tends to overestimate the losses; however, neglecting the ∼3 cm gyroradius underestimates losses.) The steady state beam ion  population is inferred from the combination of all 42 recorded beam ion populations, which is equivalent to integrating the confined distribution function in time over the run length, and implies that beam fuelling is included in the model. Steady state distributions derived using this method in cases with and without a perturbation are shown in figures 8 and 9.
A simple analytical expression is used to model the tearing mode in the simulations. The tearing mode is assumed to have created a (2,1) island chain at the q = 2 surface. In the Hamiltonian formulation of guiding centre motion, a magnetic perturbation δB is introduced through a vector potential, δB = ∇ × αB, where B is the equilibrium field and α(ψ p , θ, ζ ) is a scalar function of position. (The code uses the poloidal flux ψ p as its radial variable, with θ and ζ as the poloidal and toroidal angles, respectively.) This form is sufficient to describe the ∇ψ p structure of static magnetic perturbations. A perturbation of the form α(θ, ζ ) = α 0 sin(nζ − mθ ) results in magnetic islands of width with q = m/n and q = dq/dψ p , representing the local shear [32]. This functional form describes the tearing mode approximately [28] but is unsuitable numerically. Consider instead [5] α where x = r/a and r is the minor radius along the outer midplane (r ∼ 2ψ p /B). The quantity x 0 is the minor radius at which the mode peaks, and the quantity p ≡ m(x −1 0 − 1) is specified such that the mode peaks near q = m/n. (The value of p is rounded out to the nearest integer in order to speed up the calculation.) This form of the perturbation is a good model for the saturated, single helicity tearing mode: it extends globally  for low values of m and n; it peaks at the q = m/n surface; it vanishes at the plasma edge (i.e. conducting wall); and it has the appropriate x m dependence as x → 0. The tearing mode is a global mode and the measured island is large, so small 'bumps' in the radial eigenmode structure are not expected. The radial portion of the perturbation used in the simulations is shown in figure 5. Only a single helicity mode is modelled in the simulations. The toroidal Mirnov array clearly shows that the dominant mode is n = 1, and no higher frequency modes associated with beam driven instabilities are present. The poloidal mode number m = 2 is inferred only from the location of the phase change in the ECE data ( figure 3). Since the discharge is an elliptical, toroidal plasma, other harmonics must be present. Additional helicities could significantly increase radial transport [5,7], and neoclassical tearing modes can have higher poloidal mode numbers [18].
Determining the island width for a given perturbation strength is tricky. Though the estimate in equation (2) is not derived from the perturbation described by equation (3), it remains a good approximation when the islands are relatively small. So in the simulations, the island width for a given perturbation strength is obtained from equation (2), with α evaluated at the q = m/n = 2/1 flux surface. The island width is translated from the poloidal flux co-ordinate used by ORBIT, ψ p , to the toroidal flux co-ordinate ρ used to analyse and display data at DIII-D, and reported here as δρ. The island width can also be obtained directly by considering a Poincaré puncture plot of low energy ions (figure 10), though discerning the outline of the island is difficult. Poincaré plots compare well with estimates from equation (2), even at the largest island sizes used in the simulations, so they are reported in the results. (When the island half-width is comparable to the radial distance from the magnetic axis to the q = 2 surface, the asymmetry in α (figure 5) becomes evident as the inner radius essentially becomes fixed and only the outer radius grows with the amplitude.) The perturbation strength that gives the experimental island size δρ 0.26 is δB r /B 10 −3 , which is smaller than the experimental value calculated from the Mirnov data, δB r /B = (4.1 ± 1.6) × 10 −3 .
The current drive and reaction probabilities are inferred from the confined ion population using expressions given below. Experimental profiles of the electron and deuterium density (n e , n D ), the electron and ion temperature (T e , T D ) and the carbon density n C are utilized in the calculations ( figure 1).
Instantaneously, each beam ion drives a current proportional to where e is the charge constant, ν is the particle velocity, λ is the particle pitch, Z eff is the average ion charge and ε is the aspect ratio. The term in the curly brackets represents the effect of including electrons; the trapping factor F accounts for finite aspect ratio. In the collisionless regime, a good approximation for this factor is F (1.55 + 0.85/Z eff ) √ ε − (0.20 + 1.55/Z eff )ε [33]. The total current is the sum of the current contribution from all confined particles.
The neutron rate is based upon computations of the reaction probability per particle. From TRANSP, the neutron signal is dominated by beam-plasma reactions but ∼30% of the neutron emission comes from beam-beam reactions, with the thermonuclear contribution being negligible. For beam-plasma reactions, the reaction rate per beam ion is n D σ ν b−p , where σ ν b−p is the average value of the D-D reactivity and can be approximated by [24] Here n D and T D are the local deuterium density and temperature, σ ν is the D-D fusion reactivity [34] evaluated at the beam ion speed, and c 1 is an energy dependent coefficient that accounts for the enhancement in reactivity associated with finite ion temperature. The beam-beam reaction rate is calculated as follows. The Monte Carlo beam distribution is subdivided into spatial shells of volume V , with approximately N B 10 beam ions per shell. The reactivity σ ν b−b is calculated from the relative velocities of the ions in this radial shell, including several random values of gyrophase. The beam-beam rate from the shell is then proportional to N 2 B σ ν b−b / V . The total reaction rate due to the confined beam ions is the weighted average of the beam-plasma and beam-beam rates.

Results
The calculated loss of particles, driven current and neutron rate are shown as a function of perturbation amplitude in figure 11. To within experimental uncertainties, the calculated reductions in driven current (35%) and neutron rate (40%) are consistent with the experimental measurements when the perturbation equals the experimental island width (δρ = 0.26). Because the probability of a fusion reaction is largest near the magnetic axis and the losses are largest in the centre, the fractional particle losses in the simulations are consistently much less than the neutron rate losses. For the simulation results shown in figure 11, the particle losses scale as δB 2 r (with a correlation coefficient of 0.995 for the linear regression fit).
With the tearing mode, the distribution has about 20% fewer beam ions than without ( figure 8(a)). The tearing mode perturbation acts most effectively on passing particles, λ ≈ 1 ( figure 8(b)). The mode causes a strong reduction in central density ( figure 9). For the case with a perturbation, the steady slope of the energy distribution (figure 8(a)) shows that, though the mode is most effective at removing particles near the birth energy, it continues to cause losses as the injected particles slow down. The spatial distribution in figure 9 indicates a significant loss of fast ions at the peak location and some redistribution of particles about the peak.
The lost particles at the radial peak of the distribution cause the largest reduction in both the current drive and the neutron rate ( figure 12). The redistribution of particles intensifies the losses in both quantities. The beam-beam reactivity is most adversely affected, as almost ∼50% of beam-beam produced neutrons are lost in the simulations. Note that the largest reductions occur for ρ < 0.3, i.e. near the tearing mode.
Simulations were also performed using the classical slowing down distribution function calculated by TRANSP [16] as the initial distribution. No collisions were included in those cases. Particle trajectories were followed for ∼18 ms, and cases with and without a perturbation were considered. The calculated losses of particles, current drive and neutron rate were similar to the results of the more rigorous simulations shown in figure 11.  Figure 11. Total lost driven current (squares), neutron rate (circles) and particles (crosses) for a range of island widths. The maximum δB r /B on this plot is ∼10 −3 . Also pictured is the range in the measured island width and the neutron loss data (shaded area).

Discussion and conclusions
There are few known non-resonant transport loss mechanisms for beam ions on DIII-D, and none can really explain the data and the simulation results as well as orbit stochasticity. The helical perturbation from the tearing mode might be causing the fast ions to cross a loss boundary; however, this mechanism is not expected to be important since the beam ions are mostly passing and thus well confined. The MHD mode could also be causing magnetic field lines to become ergodic. Particles with a large component of velocity could follow these field lines as they stream around the torus, resulting in a net radial displacement and possibly enhanced losses. The energy confinement times, however, for cases with and without tearing modes are comparable, indicating little or no additional anomalous electron transport. Also, the stochastic threshold for the drift orbit surfaces is expected to be much lower than for the magnetic surfaces [4]. Considering these factors, we conclude that the tearing mode is not causing the field lines to become ergodic. The experimental observation of orbit stochasticity as a viable transport mechanism for energetic ions has some important implications. This mechanism is an additional loss channel for normally well confined passing particles; the findings strongly support the speculation in other studies that this mechanism was responsible for anomalous and enhanced losses [1]. That means that non-resonant MHD activity such as the neoclassical tearing mode, normally expected to affect only the bulk species and energy confinement, can degrade the confinement of energetic ions as well if the modes are large enough. On a positive side, orbit stochasticity has been connected with a technique of removing intermediate energy alpha particles in a reactor through the external application of helical perturbations [35]. (The removal of this 'helium ash' is necessary for a self-sustaining burn.) These results give credence to exploring this transport mechanism systematically using externally applied perturbations, which would be relevant to helium ash removal.
To summarize, this work was initially undertaken to study NBCD losses, correlated with neutron losses, in the presence of large tearing mode islands. An analytic estimate of island overlap indicates that orbit stochasticity is important for beam ions in this discharge and is the likely transport mechanism for the observed losses. Simulations were conducted to study the transport numerically. By accounting for the large beam-beam contribution to the neutron rate, the total neutron losses in the simulations are found to be consistent with measured losses at island widths comparable to the experimental value. Also, the NBCD losses are limited to ρ 0.3 in the simulations, qualitatively consistent with the data. In the simulations, the helical perturbation affects high energy, circulating beam ions, and the particle losses scale as δB 2 , as expected for orbit stochasticity. This evidence strongly suggests that intrinsic orbit stochasticity causes the transport of fast ions seen in the experiment.