Self-consistent calculation of the effects of RF injection in the HHFW heating regimes on the evolution of fast ions in toroidal plasmas

A critical question for the use of ion cyclotron range of frequency (ICRF) heating in the ITER device and beyond is interaction of fast waves with energetic ion populations from neutral beam injection (NBI), fusion reactions, and minority ions accelerated by the RF waves themselves. Several experiments have demonstrated that the interaction between fast waves and fast ions can indeed be strong enough to significantly modify the NB ion population. To model the RF/fast ion interaction and the resulting fast ion distribution, a recent extension of the full wave solver TORIC v.5 that includes non-Maxwellian effects has been combined with the Monte Carlo NUBEAM code through an RF “kick” operator. In this work, we present an initial verification of the NUBEAM RF “kick” operator for high harmonic fast wave (HHFW) heating regime in NSTX plasma.


Introduction
Self-consistent modeling of the interactions between fast waves (FWs) and fast ions, introduced either from neutral beam injection (NBI) or from fusion-generated alpha particles, is important for both present-day experiments and for ITER.The fast-ion population changes the wave propagation and absorption, while the wave damping on fast-ions modifies their distribution.The latter implies that fast-wave heating could impact and perhaps give leverage over Alfvénic instabilities driven by beam ions.Specific to NSTX [1], FWs can modify the behavior of -and, under certain conditions, even suppress -fast ion driven instabilities such as toroidal Alfvén eigenmodes (TAEs), global Alfvén eigenmodes (GAEs) and fishbones [2,3].Moreover, because of the lower toroidal field of the spherical tokamak, FW heating may accelerate fast ions to loss orbits, and this power-loss mechanism must be studied and then minimized.This paper is structured as follows: the main equations of the RF-kick operator implemented in NUBEAM are presented in Section 2. In Section 3, a brief description of the recent extension of the full wave code TORIC [4,5] is presented.In Section 4, first and preliminary tests/results are shown and discussed for NSTX plasmas.

RF-kick operator implemented in NUBEAM
The original implementation of the RF-kick operator in NUBEAM [6,7] was presented in [8,9].However, a clear verification of such implementation is still missing.In this section we briefly show the main equations of the RF-kick operator implemented in NUBEAM.The RF-kick operator e-mail: nbertell@pppl.govimplemented in NUBEAM is based on the original works by C. F. Kennel, F. Engelmann [10], T. H. Stix [11], and X. Q. Xu and M. N. Rosenbluth [12].The strong damping of FWs on fast ions occur at the ion-cyclotron resonance in the plasma and the RF resonance condition for a RF angular frequency ω rf is given by where B is the magnetic field, n is the harmonic number, k is the parallel component (to the magnetic field) of the wave-vector, E ≡ m s v 2 s /2 is the kinetic energy of the particle, µ ≡ m s v 2 ⊥,s /(2B) is the particle magnetic moment.Finally, q s and m s are the particle charge and mass, respectively.Applying a similar derivation of the Coulomb collision operator with a Monte-Carlo scheme in [12] to our problem of interest, one can evaluate the new parallel v and perpendicular component v ⊥ of particle velocity from the old velocities (v where ν rf ⊥ , ν rf , ν rf ⊥ indicate the perpendicular and parallel velocity RF diffusion processes and R 1 , R 2 are two EPJ Web of Conferences 157, 03004 (2017) independent random variables that are distributed between 0 and 1.Finally, ∆t represents the integration time step (during the particle orbiting).Moreover, the leading term of the RF quasi-linear diffusion operator operator D rf can be written as [11,13] (7) where E + is the left-handed wave electric field and J n represents the Bessel functions with argument k ⊥ v ⊥ /ω ci (k ⊥ and ω ci are the perpendicular component of the wave vector and the ion cyclotron angular frequency, respectively).In Eq. 7 and in the TRANSP implementation the righthanded wave field, E − , and the parallel wave field, E have been currently neglected.

TORIC non-Maxwellian extension
The version of TORIC implemented in TRANSP [14], which is able to solve the Maxwellian equations in the high harmonic frequency heating regimes [4], has been recently extended to include non-Maxwellian ion effects [5].The extension of the code consists mainly in the implementation of the full hot plasma dielectric tensors for arbitrary distribution function instead of the Maxwellian distribution function as originally implemented.The HHFW version of TORIC is based on the so called "Quasi-local" approximation (the reader is referred to [4] for more details).The structure of the differential form of the wave equation implemented in TORIC in the finite-Larmor-radius (FLR) approximation is maintained and 0th-order FLR coefficients of the wave equation are replaced by the corresponding elements of the full hot-plasma dielectric tensor in which the k ⊥ value in the argument of the Bessel functions is obtained by solving the local dispersion relation for FWs only.This code extension can now deal with analytical functional forms of the distribution functions (such as, anisotropic Maxwellian and slowing down distributions) and numerical distributions.In fact, the main goal of this extension was to enable the code to deal with numerical distribution provided by Fokker-Planck (for instance, CQL3D code [13]) and/or Monte Carlo particle codes (for instance, NUBEAM [6,7]).More specifically, in this work we want to use the distribution function obtained from the Monte Carlo NUBEAM code for an NSTX [1] discharge.To do that, we make use of the P2F code [15,16], which converts a discrete particle list to a 4-D continuum distribution function.The 4 dimensions are 2 cylindrical in space (R, z) and 2 cylindrical in velocity space (v , v ⊥ ).This conversion is necessary because the full hot dielectric dielectric tensor for arbitrary distribution functions, implemented in the code, takes velocity space derivatives of the distribution function, therefore function must be smooth enough for the derivatives to be robust.Several validations and tests with analytical and numerical distribution of these code have been presented in [5] and we are referring the reader to that reference.However, in Fig. 1 we show an example of the beam distribution function obtained from the P2F code starting with a NUBEAM particle list.

Preliminary test of RF kick operator
In the TRANSP code, the TORIC and NUBEAM code are currently coupled in the following way: TORIC provides the wave electric field, the perpendicular wave vector and the toroidal mode number to NUBEAM in order to run the RFkick operator and in turn, NUBEAM provides to TORIC the profiles of the flux surface average of the total energy density and the density of the beam ions.These two quantities are used to get an effective temperature given by [17] where E and n fi are the total energy density profile and the density of the beams ions, respectively.Furthermore, currently in TRANSP there is a re-normalization factor introduced to enforce the total power absorbed by fast ions evaluated by NUBEAM to match the absorbed power evaluated by the TORIC code in order to have consistency between the NUBEAM and TORIC calculations [9].This normalization factor varies on a case-by-case basis.In this paper, we present an initial verification of the NUBEAM RF "kick" operator for high harmonic fast wave (HHFW) heating regime in NSTX plasma.In particular, Figure 2 shows the NSTX discharge # 128733 used in these numerical simulations.NSTX operates with a toroidal field of 0.5 T, with typical density (3 − 10) × 10 19 m −3 , temperature T e ≈ T i < 1.5 keV.For this discharge there are two main heating systems: neutral beam injection (NBI) and HHFW.During the flattop phase 1 MW of HHFW and NBI are injected.The NB injection energy is 65 keV (see also Figure 1).The antenna phase is 90 degree, which corresponds to a toroidal wave vector k φ = 8 m −1 and the antenna frequency is 30 MHz.Full wave simulations for NSTX plasmas commonly found a significant amount of RF power absorbed by the fast ion population.In fact, several IC cyclotron resonances, where the strong interaction between FW and particles occurs, are present in the plasma, as shown in Figure 3 for NSTX discharge 128733.
The black curve in Figure 4 shows the measured value of the neutron rate for the NSTX discharge 128733,  whereas the red curve is the predicted neutron rate obtained from TRANSP without the application of the RF kick-operator in NUBEAM.One can clearly see that the level of neutron emission rate rises significantly during the HHFW heating (as indicated by the green horizontal line for reference).Because the neutron rate is dominated by beam-plasma reactions, this indicates that a significant fraction of high-energy fast ions are accelerated by the RF heating.In fact, this is commonly found in the RF simulations for NSTX plasmas when a NBI species is considered [17].If this effect is neglected, TRANSP completely fails in predicting the rise in the neutron rate.As a first test, we repeat the same TRANSP simulation includ- ing the RF kick operator and the contribution from TORIC.
The results are shown in Figure 5.We show two cases: green and red curves represent two TRANSP runs with and without the re-normalization factor between TORIC and NUBEAM, respectively.The predicted neutron rate without re-normalization factor tends to underestimate the measured neutron rate except for the very beginning of the HHFW application.On the other hand, the predicted neutron rate assuming a re-normalization factor tends to be closer to the measured one.However, it overestimates the neutron rate for about the first half of the HHFW application.In the analysis above, TORIC is using a Maxwellian distribution via the effective temperature (cf.Eq. 8) although NUBEAM provides a general fast ion distribution function.Regarding this issue, as a starting point, we have compared the power deposition profiles obtained by TORIC and NUBEAM including the RF-kick operator for an NSTX plasma described in [5].More specifically, we have assumed a Maxwellian distribution function for the fast ion population at the deposition in NUBEAM.Figure 6(a) shows the HHFW power profiles from TORIC (red and green curves) and NUBEAM (blue curve).The green and red curves are obtained assuming the fast ion density at the deposition and the fast ion density after following the particles in a single NUBEAM time step.In both cases, an effective temperature is used as shown in Eq. 8. Figure 6(b) shows the corresponding flux surface average of the fast ion density at the deposition (red curve) and after following the particles in a single NUBEAM time step (blue curve).From 6(a), the agreement between TORIC and NUBEAM is quite acceptable.The peak locations and the profiles are similar.Moreover, TORIC's result tends to be closer to NUBEAM's result when using the fast ion density after following the particles for a first time step in NUBEAM.This preliminary exercise seems to suggest that when the fast ion distribution function between TORIC and NUBEAM are similar there is a reasonable agreement between the power profiles from the two codes.However, it is important to stress that this is still a preliminary and a single result.Several other tests are underway and they will be shown in a future work.In particular, the consideration of a realistic distribution obtained starting from NUBEAM and affected by RF, will be part of a future work together with additional tests on the "RF-kick" operator in NUBEAM.

Figure 1 .
Figure 1.Fast ion distribution function at R = 1 m and z = 0 m (approximately corresponding to the magnetic axis) obtained from the P2F code starting from a NUBEAM particles list for an NSTX discharge 128733 at time = 0.199 s (pure NBI distribution without the RF effects).

Figure 4 .Figure 5 .
Figure 4. Neutron rate for NSTX discharge 128733.The black curve indicates the measured values, whereas the red curve is the predicted neutron rate obtained from TRANSP without the RF kick-operator.The green horizontal line shows for reference when HHFW is applied.

Figure 6 .
Figure 6.(a) ICRH power absorption profiles of the beam ions as a function of the square root of the normalized poloidal flux, ρ pol evaluated by NUBEAM (blue curve) with "RF-kick" operator and TORIC with the fast ion density at the deposition (red curve) and after following the particles (green curve) including fast particle losses in NUBEAM (see blue curve in Fig. b.).(b) Fast ion density at the deposition (red curve) and after following the particles including fast particle losses (blue curve) as a function of the square root of the normalized poloidal flux, ρ pol .