Monte Carlo orbit/full wave simulation of ion cyclotron resonance frequency wave damping on resonant ions in tokamaks

To investigate the experimentally observed interaction between beam ion species and fast Alfvén wave s FW d , a Monte Carlo code, ORBIT-RF f V. S. Chan, S. C. Chiu, and Y. A. Omelchenko, Phys. Plasmas 9 , 501 s 2002 dg , which solves the time-dependent Hamiltonian guiding center drift equations, has been upgraded to incorporate a steady-state neutral beam ion slowing-down distribution, a quasilinear high harmonic radio frequency diffusion operator and the wave ﬁelds from the two-dimensional ion cyclotron resonance frequency full wave code s TORIC4 d f M. Brambilla, Plasma Phys. Controlled Fusion 41 , 1 s 1999 dg . Comparison of ORBIT-RF simulation of power absorption with ﬁxed amplitudes of FW ﬁelds from TORIC4 power absorption calculation, which assumes Maxwellian plasma distributions, attains agreement within a factor of two. The experimentally measured enhanced neutron rate is reproduced to within 30% from ORBIT-RF simulation using a single dominant toroidal and poloidal wave number. © 2005 American Institute of Physics . f DOI: 10.1063/1.1935387 g


I. INTRODUCTION
In burning plasma physics experiments such as International Thermonuclear Experimental Reactor ͑ITER͒, 1 both neutral beam injection ͑NBI͒ and fast Alfvén wave ͑FW͒ in the ion cyclotron resonance frequencies ͑ICRFs͒ will be employed for auxiliary heating and current drive. In consequence, a large population of energetic ions exists in the form of injected neutral beam ions and fusion born alphas. Since energetic ions have k Ќ i ജ 1 ͑k Ќ is the perpendicular wave number and i the ion Larmor radius͒, strong damping of the FW on energetic ions even at high ion cyclotron harmonics may occur near the resonant locations satisfying rf − k ʈ ʈ = n⍀ i where rf is the wave frequency, k ʈ the parallel component of wave number, ʈ the parallel ion velocity, n the harmonic number and ⍀ i the ion cyclotron frequency. Therefore, significant high velocity ion tails can be produced in the plasma when the FW is launched into the plasma. The presence of energetic ion species can modify the plasma behavior in important ways such as magnetohydrodynamic ͑MHD͒ stabilization or destabilization and momentum transport/plasma rotation by the ICRF generated energetic tails. These issues relating the interaction between the ICRF wave and energetic ions must be well understood for successes of future tokamak experiments and achievement of ITER physics goals.
Answers to some of the issues might come from understanding previous experiments in existing tokamaks. Beam ion absorption of the FW at moderately high cyclotron harmonics has been observed experimentally in JT-60U ͑Japan Tokamak-60 Upgrade͒, 2 JET ͑Joint European Torus͒, 3 TEXTOR ͑Tokamak Experiment for Technology Oriented Research͒, 4 and DIII-D. [5][6][7][8][9] In the DIII-D tokamak, absorption of the FW on beam ion species has been observed during fast wave current drive ͑FWCD͒ experiments in deuterium beam injected plasma using either two simultaneous frequencies ͑60 and 83 MHz͒ or separate frequencies from 60 MHz to 117 MHz. These DIII-D experiments have demonstrated that the presence of ion cyclotron harmonic resonances in the plasma can lead to significant damping of the FW on energetic beam ions generated during neutral beam injection ͑NBI͒ when the cyclotron harmonic layers ͑n =4 ϳ 8͒ are near the magnetic axis. Most recently, experiments in the National Spherical Torus Experiment ͑NSTX͒ have reported 10 significant interaction between high harmonic fast wave ͑HHFW͒ and energetic beam ions. As evidence, peaked pressure profile of beam ions has been observed in the central region of plasma with a significant enhancement of measured neutron rate when the HHFW was launched to NB preheated discharges. The ion energy spectrum from neutral particle analyzer ͑NPA͒ showed also strong enhancement of energetic tails on beam ion distribution above the injected beam ion energy, which leads to enhanced neutron rate.
Conventionally, full wave 11,12 or ray-tracing codes 13,14 combined with Fokker-Planck code assuming zero banana width have been used to study the interaction between the FW and the plasma. Previous numerical analysis on the experimental measurements in DIII-D tokamak has been also done using similar simple Fokker-Planck model 9 and the ICRF code Power deposition for ION cyclotron resonance heating ͑PION͒. 8 However, these approaches do not take into account both finite orbit width of non-Maxwellian ions produced from the radio frequency ͑rf͒ interaction and the radial scattering of energetic particles across flux surfaces by collisions and waves. For example, since conventional full wave codes assume local Maxwellian distributions for plasma species, the predicted power absorption of the FW on ion species may depend significantly on the initially assumed distribution. In addition, since rf excited ions may also go through large orbit excursions, absorbed power deposition profiles on ion species may be much broader or their power fractions may be reduced due to orbit loss to a wall. In a large device like ITER, this finite orbit effect is expected to be important since several resonance layers can be located close to each other spatially across the minor radius. Therefore, these effects should be included for more reliable predictions of future DIII-D experiments and burning plasma physics.
A similar study on the effects of finite orbit width and radial diffusion of rf-induced non-Maxwellian ions due to the interaction between the FW and plasma ion species is the work by Hellsten et al. 15 It used the SELFO code that calculates self-consistently the wave field from a global full wave code LION 16 and the non-Maxwellian velocity distribution of fast ions from a Monte Carlo code FIDO 17 solving the Langevin equations for the invariants of motion of the unperturbed orbits for each resonant ion species. 15 The dielectric tensor is constructed from the calculated distribution function.
ORBIT-R 18 differs from SELFO in that it follows the Hamiltonian guiding center orbit drift motion for each ion species. Main features of ORBIT-RF are as follows. It solves the Hamiltonian guiding center drift equations to follow trajectories of test ion species in 2D axisymmetric numerical magnetic equilibrium under Coulomb collisions and ICRF quasilinear heating. Monte Carlo collision operators for pitchangle scattering and drag calculate the changes of test ions in parallel velocity and pitch angle due to Coulomb collisions between test ions and background plasma. A rf-induced random walk model describing the fast ion stochastic interaction with the FW is implemented to reproduce quasilinear diffusion in velocity space, assuming that test ions lose their phase information with the FW through successive collisions before they re-enter the resonance region. 19,20 A generalized high harmonic rf diffusion operator 20 calculates perpendicular rf kicks at resonances including Doppler shifts. The two dimensional full wave code TORIC4 12 is coupled to ORBIT-RF to determine radial profiles of the FW field and its wave numbers in the plasma. TORIC4 solves the wave equations in ion cyclotron range of frequencies in 2D numerical magnetic equilibrium for a given toroidal mode number ͑n ͒ where dielectric tensors in a hot plasma are calculated using modified Bessel functions ͓I n : n is a harmonic number ͑0,1,…͔͒ and its derivatives at all orders in Larmor radius to include nth ͑n ӷ 2͒ harmonic resonance interaction in wave equations. In this paper, we did not attempt to self-consistently couple the wave fields from TORIC4 and the non-Maxwellian ion distribution from ORBIT-RF, which will be left for future study. In TORIC4, a dielectric tensor is constructed using Maxwellian distribution of plasma species. Therefore, the magnitudes of FW fields are fixed during simulation time. Steady-state slowing-down distribution of beam ion species is modeled using a re-injection method of thermalized beam ions. These new features of ORBIT-RF will be described in detail in Sec. III. Even though several harmonic cyclotron resonance layers of resonant ions exist simultaneously in the plasma in actual experiments, they are sufficiently well separated to justify modeling only a single cyclotron resonance layer in this work. The extensive physics features of ORBIT-RF are capable of identifying the experimentally demonstrated strong interaction between the FW and injected beam ions in the DIII-D tokamak. The simulation results show that the physical insight and predictions provided by this approach qualitatively explain the experimental observations and point to further improvements needed in the models to simulate experiments more quantitatively.
This paper is organized as follows. The experimental conditions on DIII-D for high harmonic ion cyclotron damping of the FW on injected beam ion species are reviewed in Sec. II. Section III describes the main features of ORBIT-RF: modeling of steady-state slowing-down distribution of beam ions, analytical expressions for quasilinear rf diffusion operator with arbitrary harmonic absorption, Coulomb collision operators between test ions and background plasma, and numerical expressions for neutron enhancement and absorbed power on resonant ion species as numerical diagnostics. In Sec. IV, the numerical results are discussed for a selected DIII-D discharge. Here, detailed modeling from the full wave code TORIC4 is first presented for a numerical equilibrium corresponding to the chosen DIII-D discharge together with results from sensitivity study of TORIC4. In sequence, the numerical results from ORBIT-RF coupled with wave fields from TORIC4 are compared with power absorption calculated with TORIC4 for Maxwellian distributions. Lastly, ORBIT-RF predictions are compared with the DIII-D experimental measurement in terms of absorbed powers for resonant ion species and neutron enhancement. A summary and plans for future work are presented in Sec. V.

II. EXPERIMENTAL CONDITIONS IN THE DIII-D TOKAMAK
Among many experimental results in the DIII-D tokamak [5][6][7][8][9] using different frequency ranges ͑60-117 MHz͒ at different harmonics ͑fourth-eighth͒, we focus on the experiments performed during 1998 campaign for fourth harmonic central damping of 60 MHz FW on injected deuterium beam ions. 9 The major radius was R 0 = 1.7 m, the minor radius a = 0.6 m, the usual toroidal magnetic field B 0 = 1.9 T and the plasma current I p = 1.2 MA. Four-element phased array were used to launch the 60 MHz FW with a peak in the vacuum spectrum at n ʈ = k ʈ c / Ϸ 5. Deuterium beam ions of 70-86 keV are injected in the co-current direction with P NB = 1.6-2.7 MW, which provides a fast ion population in the plasma. After the NB preheated plasma reached a steady state, rf power with P rf = 1.0-1.5 MW was coupled for a 1 -2 s pulse. The central electron density n e ͑0͒ was typically 2.0-5.0ϫ 10 13 cm −3 , the electron temperature T e ͑0͒, 2.0-4.0 keV, the plasma ion temperature T i ͑0͒, 2.0-3.0 keV, and the effective charge Z ef f Ϸ 2.0. At B 0 = 1.9 T, the 60 MHz FW interacts with deuterium beam ions at three resonance locations, as shown in Fig. 1. The fourth harmonic resonance layer is located near plasma center.
Among three vertical charge exchange lines of sight, V1, V2 and V3, in Fig. 1, only V2 located near the fourth harmonic deuterium resonance layer detected significant energetic tails during rf pulse, 9 which confirms the formation of energetic beam ion tail population in the central region. As a consequence, a peaked ion pressure profile due to the increase of perpendicularly stored energy was observed experimentally with enhanced d-d neutron emission in the central region of plasma on beam ion slowing-down timescales. The measured neutron rate enhancement ͓Eq. ͑10͔͒ was mostly between 1 and 2 during this campaign. This observation strongly suggests significant damping of the FW on energetic beam ion species at the fourth harmonic layer. The experimental observations described in this section will be investigated through numerical simulations using ORBIT-RF coupled with a full wave solver, TORIC4, in Sec. IV.

A. Modeling of steady-state slowing-down distribution of beam ions
The injected neutral particles are assumed to be all ionized by background plasma. Initially ionized beam ions are distributed randomly according to n f ͑͒ ϰ ͑1− 2 ͒ 2 in real space while they have an initially monoenergetic distribution in velocity space where n f and are a beam ion density and normalized poloidal flux, respectively. In this work, only the full energy component is considered. Pitch angle of beam ion is calculated by = ʈ / = ͑B T R T ͒ / ͑BR͒ ͑Ref. 21͒ where ʈ is the parallel velocity of beam ion, the velocity of beam ion and R T the tangency radius. To represent real numbers of injected beam ions at any time consistent with the injected beam power, a time-dependent weight is assigned to each test beam ion. For simplicity, we assume that test beam ions initially have the same weights, which are determined by the beam ion source rate S calculated using the given P NB ͑NB injection power͒ and E b ͑beam ion energy͒. Therefore, each beam ion has a distribution at each time in ͑E b , , , , w͒ space where and w are the poloidal angle and weighting, respectively.
Injected monoenergetic beam ions lose their energy and momentum to background plasma ions and electrons through pitch-angle scattering and Coulomb collisions, and eventually are thermalized. To model a steady-state slowing-down distribution of injected beam ions, we re-inject thermalized beam ions periodically into the plasma. Reinjection times are decided initially to be several ms consistent with thermalization time. While beam ions are slowing down, they are monitored at every reinjection time whether their energies have dropped below a critical value ͑typically chosen to be 1.5 T i where T i is the temperature of thermal plasma ions͒. If so, new beam ions are injected at their birth energy with a source rate compatible with P NB . Then, the weighting of reinjected beam ion is readjusted to be consistent with constant input power, and simultaneously with ion fueling rate of beam. This procedure is carried through a few slowingdown times of beam ions. Results show that this approach, though simple, closely reproduces a steady-state slowingdown distribution of beam ions within a few slowing-down times.

B. Quasilinear high harmonic RF diffusion operator
When ions pass through ion cyclotron resonance layer satisfying the condition rf − k ʈ ʈ = n⍀ i ͑n: harmonic num-ber͒, they may either absorb or lose energy from the wave depending on their phase with polarization of the FW. Since ion phase information is lost through successive collisions 19 before they re-enter the resonance region, the phase between ions and the wave can be assumed to be random. Therefore, we introduce a simplified rf-induced random walk model to reproduce quasilinear diffusion in velocity space through stochastic interaction with the FW. The interaction of ions with the FW at Doppler shifted resonances changes the velocity of ions in both parallel and perpendicular directions. In this work, the change of ʈ due to finite parallel wave number ͑k ʈ 0͒ is ignored since a significant change is made mostly in perpendicular component. Therefore, the increment of magnetic moment ͑⌬ rf ͒ ͓Eq. ͑1͔͒ due to the rf interaction is obtained by calculating a mean value ͓Eq. ͑2͔͒ and random variance ͓Eq. ͑3͔͒ from quasilinear equation governing rfinduced particle diffusion in velocity space: R s is a random number between 0 and 1 where the factor 2 ͱ 3 is such that where Ai is the Airy function in order to include the case when the successive interactions between a particle orbit and the cyclotron resonances are close to each other. In Eq. ͑2͒, ͉E + ͉ is the absolute magnitude of the left-hand polarized component of wave electric field rotating in direction of ions. The contribution from right-hand component ͉E − ͉ is ignored in this work. To determine the magnitudes of ͉E + ͉, k Ќ and k ʈ in Eqs. ͑2͒ and ͑3͒, the 2D full wave code TORIC4 is coupled into ORBIT-RF. The radial profiles from TORIC4 are passed on to ORBIT-RF to calculate the increase of perpendicular energy of ions passing through resonance surface. The detailed approach is described in Sec. IV.

C. Collision operators
The change in velocity of test ions due to Coulomb collisions with plasma ions and electrons is calculated at each time interval ⌬t using the slowing-down frequency between ions-ions and ions-electrons, 23 The change in pitch angle induced by scattering between test ions and plasma ions is modeled by 23 is the critical energy, A f the atomic mass number of test ion, ⌳ f is the Coulomb logarithm, E f is the test ion energy, n i is the plasma ion density, n e is the plasma electron density and the subscript f denotes test ions and i for background thermal ions. We assume that the density and temperature of background plasma ions and electrons are fixed during simulations. Their initial profiles are reproduced from experimental data at a given time.
D. Neutron enhancement "S n … and absorbed power "P abs … As a convenient measure of the FW damping on beam ions, neutron enhancement S n is calculated using the ratio of neutron reaction rates between NB only ͑no ICRF͒ and NB coupled with the FW 9 given by where n c is the number of Monte Carlo test ions, ͗͘ the reaction rate for beam plasma and w j the weighting of each test ion.
Power absorption on resonant ion species due to rf kicks is calculated by where ⌬E Ќ is the change of perpendicular energy due to rf kicks at each time step, and t run is simulation time of the ORBIT-RF ͑typically several slowing-down times͒. The power balance between P abs using Eq. ͑11͒ and power transferred by energetic ions to background is cross checked.

IV. ORBIT-RF NUMERICAL RESULTS
The DIII-D discharge 96043 at 1920 ms is selected to investigate the experimental observations described in Sec. II. The background plasma consists of deuterium majority ions, hydrogen minority ions and electrons. Corresponding radial profiles for density ͑n e ͒ and temperatures ͑T e and T i ͒ are shown in Fig. 2. In Table I 96043.1920 are described. With f rf = 60 MHz at B axis = 1.9 T, the fourth harmonic resonance layer of deuterium beam ions ͑equivalent to second harmonic of H minority ions͒ is located at = 0.017 ͑R ϳ 164 cm͒, which is several cm inboard from the magnetic axis, as shown in Fig. 1.

A. Modeling of DIII-D experiment with TORIC4
1. Approach used for TORIC4 modeling Figure 3͑a͒ shows an antenna toroidal power spectrum as a function of n ʈ in vacuum for a typical DIII-D 60 MHz FWCD experiment. The integrated power in the range of 3.5ഛ n ʈ ഛ 5.7 ͑equivalent to 11ഛ n ഛ 17: n is the toroidal mode number͒ is approximately 40% of total integrated power from all n ʈ contributions in full toroidal spectrum. In TORIC4, the antenna is modeled by a current sheath with given antenna length l A on magnetic surface = A . For the DIII-D antenna, l A is set to be ϳ44 cm. Resulting antenna current spectrum as a function of poloidal Fourier mode number ͑m͒ is displayed in Fig. 3͑b͒, showing a symmetric spectrum with a peak at m = 0. Even though actual poloidal spectrum for the DIII-D antenna 60 MHz is rather peaked ͑not symmetric͒; 7 in this paper, we run TORIC4 using this simplified symmetric poloidal spectrum due to the lack of accurate antenna model in TORIC4. Based on the relationships k 2 = k r 2 + k 2 + k ʈ 2 = 2 / A 2 and ͑m / r͒ 2 ͑=k 2 ͒ ഛ 2 / A 2 where A is Alfvén velocity, we can estimate approximate number of poloidal mode number ͑N m ͒ for the relatively short antenna length ͑a half length is ϳ22 cm͒. With parameters shown in Fig. 2 and Table I, ͉N m ͉ is calculated as roughly 20. There-fore, the poloidal mode number that we have selected in this paper ͑N m = 31: −15ഛ m ഛ 15͒ seems reasonable, assuming negligible contributions from high poloidal mode numbers ͉͑m͉ Ͼ 15͒.
To identify well-converged wave solutions from TORIC4, robustness and sensitivity of the code TORIC4 to some input parameters have been tested out extensively. We found that TORIC4 demonstrated good convergences in general. Figure 4 shows that total absorbed powers ͓P ͑MW/ kA 2 ͔͒ by plasma species as a function of radial grid mesh ͑N r ͒ and radial profiles of real E + ͑V/m͒ along equatorial plane for different N r at a given n = 11 for two different plasmas.  Fig. 5͑a͒. It is clearly seen from Figs. 4͑b͒ and 4͑d͒ that the existence of beam ions in the plasma changes the radial structure and magnitude of wave field. The differences in total absorbed power and magnitude of wave field from different N r are less than a few %. Therefore, all TORIC4 simulations presented in this paper have been done with the resolution of 200N r ϫ 31N m . Approximately   seven wavelengths are shown in Figs. 4͑b͒ and 4͑d͒. This is confirmed by simple calculation using DIII-D parameters where a wavelength is calculated to be about 9 cm when a plasma radius is ϳ62 cm.

TORIC4 simulation results for a specific DIII-D discharge
To identify the effect of existing beam ions on parasitic ion power absorption from the full wave code TORIC4, we run TORIC4 for two cases with given rf input parameters in Table  I, numerical equilibrium for Shot No. 96043.01920 and plasma kinetic profiles shown in Fig. 2. The first case ͑Case 1͒ is for a plasma without beam ions ͓thermal D ͑95%͒minority H ͑5%͔͒ where only second harmonic resonance heating on thermal H minority ions is available. The second case ͑Case 2͒ is when the plasma consists of thermal D ͑88%͒-minority H ͑5%͒-Beam D ͑7%͒ where both second harmonic heating on H minority ions and fourth harmonic heating on D beam ions are included. To estimate averaged power absorption by resonant ion species over the dominant toroidal mode spectrum for each case, we assume same power weighting for each n in the range of 11ഛ n ഛ 17, instead of taking into account more realistic weightings of power for different n shown in Fig. 3͑a͒. Table II shows power absorptions by H minority ions at second harmonic resonance and electrons for each n for Case 1. Very negligible damping on background D ions is predicted while ϳ16-39% of the ICRF power is damped on electrons. The TORIC4 simulations predict averaged power absorbed by the H minority of 0.33 MW ͑ϳ75% of ICRF power͒ for the range of 11ഛ n ഛ 17.
In Table III, predicted power absorptions by three different ion species for Case 2 when beam ions exist in the plasma are summarized for a given n = 11. Here, three different local Maxwellian temperature profiles are assumed for beam ions as shown in Fig. 5. When beam ion species exist, negligible dampings on electrons and majority background ions are found. Beam ion absorbed power predicted from TORIC4 is in the range 30-75%, depending on initial local Maxwellian beam ion temperature. We find that power absorption on beam ions increases as beam ion temperature is getting higher.

Approach for coupling of TORIC4 full wave solution to ORBIT-RF
In order to couple the TORIC4 results to the ORBIT-RF, we introduce further approximations to simplify simulations. First, the TORIC4 simulations presented in previous subsection store wave solutions for each Fourier poloidal mode number ͑m: total 31 Fourier poloidal modes͒ along an equatorial plane for a given n . Instead of doing simulations with each wave field component of different m, we do simulations only with the biggest wave field component from dominant poloidal mode number for a given n . Therefore, we assume that rf power ͑ϳ0.44 MW for a given n ͒ is transferred by a single toroidal and poloidal mode.
Second, the radial magnitudes of ͉E + ͉ from TORIC4 are rescaled using a given rf power since TORIC4 calculates the magnitudes of ͉E + ͉ for 1 A antenna current. Our approach is to recalculate ͉E + ͉ to get realistic numbers corresponding to total input power, using the ratio of total power absorption from experiment ͑P EXP ͒ and that from TORIC4 ͑P TORIC ͒. For example, the TORIC4 simulations resulted in total absorbed power on plasma species, P TORIC = 3.19ϫ 10 −1 MW/ kA 2 for Case 1 ͑no beam ion͒ and 9.61ϫ 10 −1 MW/ kA 2 for Case 2 ͑with beam ion͒ for a given n = 11. Assuming P EXP = 0.44 MW for n = 11, we calculate current I ant in antenna as I ant ϳ ͱ 0.44/ 3.19ϫ 10 −1 ϳ 1.2 kA for Case 1 and I ant ϳ ͱ 0.44/ 9.61ϫ 10 −1 ϳ 0.7 kA for Case 2. The magnitudes of ͉E + ͉ from TORIC4 is multiplied by I ant to obtain the wave field amplitude for ORBIT-RF simulation. Figures 6 and 7 show radial profiles of rescaled magnitudes of ͉E + ͉ ͑V/m͒ along equatorial plane at the resolution of 200N r ϫ 31N m for a single dominant wave mode for Case 1 ͓Fig. 6͑a͔͒ and Case 2 ͓Fig. 7͑a͔͒. Corresponding radial profiles of n ʈ and n Ќ 2 for each case are also displayed in Figs. 6͑b͒ and 6͑c͒ and Figs. 7͑b͒ and 7͑c͒, respectively. These renormalized radial profiles are passed on to the ORBIT-RF to calculate the increase of perpendicular energy of resonant ions passing through resonance surface and thus absorbed power. Figure 8 shows ͉E + ͉ ͑V/m͒ as a function of ͑X,Z͒ for Case 2 before it is rescaled with estimated total antenna current. Coupling of ͉E + ͉ as shown in Fig. 8 with the ORBIT-RF has not been attempted yet in this work. Instead, the two-dimensional refractive effect of incident wave field in the plasma is implemented using a simple wave cone model inside the plasma in summary.

The ORBIT-RF simulation results
We follow the trajectories of test ions by solving Harmiltonian guiding center drift equation at each time step ͑⌬t is usually ϳ10 −5 s͒ using a Runge-Kutta fourth order integra- tion scheme subject to the FW heating and slowing-down collisions/pitch angle scattering with background ions and electrons. Simulation time is typically several slowing-down times of test ion species, which requires usually a few ten thousand transit times of test ions. Grid size used is ͑ , ͒ = ͑100, 50͒. In ORBIT-RF, electrons are considered as only background species neutralizing charge. For simulations presented in this paper, we used mostly 10 000 particles to rep-resent Monte Carlo test ions. Simulations had been performed on LINUX cluster FI at GA using 20 processors. First, two ion species plasma ͑Case 1͒ is simulated using the ORBIT-RF coupled with the biggest electric field component from the TORIC4 wave solutions for each n mode ͑11 ഛ n ഛ 17͒ using the approach described in Sec. IV. For example, Fig. 6͑a͒ shows the TORIC4 wave solution for ͉E + ͉ for n =11/m = −1. To simulate resonant interaction between hydrogen minority ions and the FW, we set minority ions as test ion species. Initial distribution of test minority ions in velocity space is assumed to be Maxwellian with ͗T e ͘ϳ2.7 keV on magnetic axis. In real space, they are distributed with a probability n j ͑͒ ϰ ͑1− 2 ͒ 2 . The second harmonic resonance of minority ions is located at = 0.017. To model a steadystate regime, thermalized ions below a minimum energy ͑E min = 1.5 keV͒ through Coulomb collisions after excited by the FW are re-Maxwellized and then reinjected into the plasma periodically during a given simulation time. The time period for re-Maxwellization and reinjection is determined by the ratio of pitch angle scattering time between ions and transit time of thermal ions, which is approximately 1000 transit times of 2.7 keV thermal ions. Figure 9 shows predicted power absorption by H minority ions at second harmonic resonance from ORBIT-RF for  each n . Absorbed power calculated using Eq. ͑11͒ ranges from ϳ0.45 MW to ϳ0.86 MW depending on n at the end of 20 000 transit time of test ions. Assuming the same weightings of power for each n to rescale wave field magnitude, we can estimate roughly averaged power of ϳ0.63 MW in the range of 11ഛ n ഛ 17 from the ORBIT-RF simulations, which is within a factor of two of the TORIC4 result of ϳ33 MW ͑Sec. IV͒. Figure 10͑a͒ shows initial Maxwellian energy distribution of minority thermal ions before rf on. Energy distribution of excited ions at the end of 20 000 transit time after rf on is shown in Fig. 10͑b͒. Energetic tails of minority ions are extended up to some hundred keV. At this time, power balance between the difference of the calculated power going to energetic ions using Eq. ͑11͒ and the power transferred by energetic ions to the background is within 0.01%. Since the ion distribution from ORBIT-RF is highly non-Maxwellian, one should not expect exact agreement in power absorption between ORBIT-RF and TORIC4.
Next, we simulate resonant interaction between injected deuterium beam ions and the FW using the ORBIT-RF coupled with the TORIC4 wave solutions for a given n = 11 as shown in Fig. 7. For this simulation, we set beam ions as test ion species. Initial distribution of beam ions in velocity space is assumed to be monoenergetic with a full energy component of 80 keV. In real space, they are distributed with a probability n j ͑͒ ϰ ͑1− 2 ͒ 2 . Detailed description about slowingdown distribution of beam ions was given in Sec. III. To model a steady-state regime, thermalized beam ions are reinjected at their birth energy ͑80 keV͒ into the plasma periodically during a given simulation time. The fourth harmonic resonance of deuterium beam ions is located at = 0.017. To model actual NB+ ICRF experimental situation realistically, NB only heating is first turned on for ϳ60 000 transit times ͑a few slowing-down times of 80 keV deuterium beam ions͒. Almost all beam ions have gone through at least one thermalization within this duration to reach a steady-state distribution. Figure 11 shows a spatial distribution of beam ion temperature ͓Fig. 11͑a͔͒ and energy distribution of beam ions ͓Fig. 11͑b͔͒ at the end of ϳ60 000 transit times right before rf on. The plot with ϩ in Fig. 11͑b͒ is calculated from an analytical expression for a steady-state slowing-down distribution of beam ions. 24 Then the FW heating is turned on for an additional ϳ30 000 transit times to explore the interaction of beam ions with the FW. The three TORIC4 cases with different Maxwellian beam ion temperature profiles discussed in the previous section ͑Fig. 5͒ have been simulated separately with ORBIT-RF by coupling each TORIC4 wave solution with ORBIT-RF. ORBIT-RF predicts that beam ions absorb ϳ0.32 MW for Case ͑a͒, ϳ0.39 MW for Case ͑b͒ and ϳ0.40 MW for Case ͑c͒. We find that the TORIC4 predictions are getting close to ORBIT-RF results as higher initial beam ion temperature is assumed in TORIC4 simulations. The large discrepancy between TORIC4 and ORBIT-RF for Case ͑a͒ are expected from the fact that the interaction between energetic beam ion tails and the FW are not being considered in TORIC4.
In Fig. 12, spatial distributions of deuterium beam ions are displayed at t =0,30 000,60 000,90 000 transit times. The FW is turned on from t = 60 000 transit times to 90 000 transit times. It is seen that the energetic tails extend up to a few hundred keV above beam injection energy ͑80 keV͒ after the ICRF heating is turned on. These tails are built up in relatively broad layers ͑0 Ͻ Ͻ 0.4͒ encompassing Doppler shifted fourth harmonic resonances as shown in Fig. 12͑d͒. In Fig. 13, radial profiles of absorbed power density are compared between the NB only and NB+ ICRF cases. Figure 14 shows that significant change is made in the perpendicular velocity, as explained in Sec. III. Recently, more quantitative energetic particle distribution was measured on NSTX 10 when high harmonic ICRF was applied to a neutral beam heated plasma, which again showed the characteristic high energy tail to several hundred keV which is again consistent with ORBIT-RF simulation.

C. Comparison of ORBIT-RF simulation results with experimental measurements
In this section, the ORBIT-RF simulation results presented in Sec. IV B are compared with experimental results in the DIII-D tokamak. Figures 15 and 16 show experimental measurements from discharge 96043 and other similar L-mode discharge plasmas, 9 demonstrating a strong interaction between beam ions and FW. Even though the neutral particle analyzer ͑NPA͒ used in the 1998 campaign was unable to measure neutrals above the beam injection energy, the increase in energetic ion population over the case with beam only is clearly shown in Fig. 15͑a͒. It is notable that the rf heating affects higher energy beam ions more strongly than lower energy ions. Figure 15͑b͒ shows radial profiles of measured thermal pressure including electrons ͑solid͒, calculated classical beam ion pressure ͑dotted͒ and adjusted beam pres-sure ͑dashed͒ that yields properly aligned ECE data. The increase of beam ion pressure with rf heating is due to the increase of perpendicular stored energy of beam ions through the interaction with rf wave. Enhanced neutron reaction rates during rf heating, as shown in Fig. 16, illustrate that deuterium beam ions are accelerated. Discharge 96043.01920 demonstrated enhancement of measured neutron rate S n ϳ 1.4.
The experimental measurements described above are qualitatively reproduced from the ORBIT-RF simulations presented in Sec. IV within a reasonable agreement. Under the assumption of no loss of energetic ions such as orbit loss and charge exchange loss, the ORBIT-RF simulation results in neutron enhancement factor of S n ϳ 1.8 ͓using Eq. ͑10͔͒, which is within 30% of the experimental value ͑S n ϳ 1.4͒. In Fig.  17, the radial profiles of beam ion pressure calculated from ORBIT-RF are compared between NB only and NB+ ICRF where the contributions from thermal electrons and ions are not included. Contrast to experimental measurements in Fig.  15͑b͒, a smaller difference between NB only and NB + ICRF is predicted which can be ascribed to the fact that the ORBIT-RF simulations presented in this paper were done by assigning only 40% of total experimental input rf power to a single toroidal/poloidal mode number. A more comprehensive simulation planned for the future which accounts for all the input power should show a larger difference.

V. SUMMARY AND FUTURE PLAN
A Monte Carlo code ͑ORBIT-RF͒ is coupled with the ICRF full wave code ͑TORIC4͒ to investigate experimentally observed interaction between injected beam ion species and the FW at high harmonics in a selected DIII-D discharge. The ORBIT-RF coupled with TORIC4 wave solutions using a single dominant toroidal and poloidal wave number demonstrates the capability to validate the qualitative consistency of simulations with experimental observations. The simulation results using a previous DIII-D FWCD discharge validated the interaction between deuterium beam ions and FW at the fourth harmonic resonance. Predicted power absorption from ORBIT-RF is benchmarked against TORIC4 calculation with a high temperature Maxwellain distribution to model the beam ions. Comparison results between the prediction from TORIC4 and that from ORBIT-RF point out that the interaction between a self-consistent energetic beam ion tails and FW may be important for a more quantitative prediction of the power absorption in experiments. Estimated enhanced neutron rate ͑ϳ1.8͒ from this specific simulation is reproduced within 30% of experimental measurement ͑ϳ1.4͒. Tail energies of injected beam ions are formed up to a few hundred keV above injected beam energy with a relatively broad spectrum when FW is turn on, which comes from the finite orbit effect of resonant beam ions circulating around resonance location near the magnetic axis and Coulomb collisions which scatter the energetic ions radially. Figure 18 shows the radial profiles of absorbed power density between TORIC4 and ORBIT-RF. A much broader power spectrum is found from the ORBIT-RF simulation which comes from orbit effect.
To consider the two-dimensional refractive effect of full wave solution, we have also modeled a wave cone in the plasma with r 0 = 0.2 ͑normalized amplitude on axis͒ and 0 = ± 30°from equatorial plane in poloidal cross section. The ions outside this wave cone are excluded from the interaction with the FW. Power absorption with this wave cone model is ϳ15% less than that of no cone model treated in this paper. Therefore, it is expected that coupling the 2D full wave solution with ORBIT-RF as shown in Fig. 8 would provide more quantitative prediction. Further self-consistent simulations should include feeding back the resulting non-Maxwellian distributions of beam and minority ions from ORBIT-RF to TORIC4 in an iterative manner to update full wave calculations. This work will help address the effect of selfconsistent interaction of wave propagation and particle orbit trajectories including Doppler shifts. Further implementation of additional energetic ion loss mechanisms is also planned in the near future.