A Finite-Orbit-Width Fokker-Planck solver for modeling of energetic particle interactions with waves, with application to Helicons in ITER

. The bounce-average (BA) finite-difference Fokker-Planck (FP) code CQL3D [1,2] now includes the essential physics to describe the RF heating of Finite-Orbit-Width (FOW) ions in tokamaks. The FP equation is reformulated in terms of Constants-Of-Motion coordinates, which we select to be particle speed, pitch angle, and major radius on the equatorial plane thus obtaining the distribution function directly at this location. Full-orbit, low collisionality neoclassical radial transport emerges from averaging the local friction and diffusion coefficients along guiding center orbits. Similarly, the BA of local quasilinear RF diffusion terms gives rise to additional radial transport. The local RF electric field components needed for the BA operator are usually obtained by a ray-tracing code, such as GENRAY, or in conjunction with full-wave codes. As a new, practical application, the CQL3D-FOW version is used for simulation of alpha-particle heating by high-harmonic waves in ITER. Coupling of high harmonic or helicon fast waves power to electrons is a promising current drive (CD) scenario for high beta plasmas. However, the efficiency of current drive can be diminished by parasitic channeling of RF power into fast ions, such as alphas, through finite Larmor-radius effects. We investigate possibilities to reduce the fast ion heating in CD scenarios.


General description of FOW features
An accurate modeling of fast ion interaction with RF waves should include Finite-Orbit-Width (FOW) effects. The bounce-average (BA) finite-difference Fokker-Planck (FP) code CQL3D [1] has been recently upgraded to incorporate all necessary FOW physics [2]. The general formulation is based on writing the original 6D FP equation in canonical action-angle space, averaging over periodic angle variables, and making a transformation from the 3D canonical action space to a computationally convenient set of invariants I = (I 1 , I 2 , I 3 ). Our choice is dictated by the desire to keep computational grids as close as possible to the original Zero-Orbit-Width (ZOW) version of CQL3D, and to maintain direct visualization of the results in velocity and configuration space. Thus, we adopt I = (u 0 ,  0 , R 0 ), where R 0 is the major radius coordinate along the minimum-B surface (the midplane of tokamak, in case of up-down symmetrical equilibrium). For each given R 0 point, the value of particle speed (relativistically, momentum-per-mass) u 0 and value of the pitch-angle  0 at this point determine a unique orbit.
The transformation from canonical action space J to our I = (u 0 ,  0 , R 0 ) space results in the general form of the bounce-averaged FPE for the particle distribution function f 0 (I,t) [2]: where   |J/I| is the Jacobian of the transformation, and the BA diffusion coefficients and advection terms are expressed through the local diffusion D uu (from collisions and/or resonant RF heating) and through the local collisional friction F u , The bounce-averaging involves the transformation coefficients I/u, which effectively give rise to nonlocal neoclassical phenomena, such as the radial diffusion described by D 33 = (R 0 /u)D uu (R 0 /u). Also, Eq. (1) includes the BA particle source operator S 0 , such as from NBI or fusion alphas. The initial benchmarking of the code is performed for calculation of the bootstrap current [2]. The results from other tests will be published elsewhere.
In this paper, we apply the FOW version of the code to simulation of alpha-particle heating by high-harmonic fast waves in ITER. Using high frequency fast waves for electron current drive (CD) is a promising scenario for high beta plasmas [3]. However, its efficiency can be diminished by parasitic channeling of RF power into fast ions, such as alphas, through wave-particle resonant Larmor-radius effects. Here, we consider two ITER scenarios [4,5], one with a relatively low content of alphas, P -source = 58 MW, and another with a high content of alphas corresponding to P -source = 112 MW.
The parameters that we change in simulations are the wave frequency and the vertical position of antenna.

ITER scenario LH_60s_3
This scenario [4], which is based on the total auxiliary power of 73 MW (33 NB + 20 EC + 20 LH) has a high central temperature T e0  T i0  31 keV and high density n e0  8.5e19 m -3 (n D = n T = 3.2e19 m -3 ), which results in a high power of fusion alphas P -source = 111.6 MW. The other plasma parameters are I p = 12.5 MA, B 0 = 5.3 T. The fast waves at 500 MHz (called "high harmonic" or "helicon" waves in this frequency range) can encounter a wide range of cyclotron harmonics for alphas, / cα = 9-17. As the first step of the modeling we use the GENRAY ray-tracing code [6] to obtain linear absorption on ions and electrons, RF field components and other data along rays, which is further imported by the CQL3D code to form the bounce-average quasilinear RF diffusion operator for the heating of alpha-particles (the Fokker-Planck'd species). As the fusion-born alphas are slowing down by collisions on one hand, and absorbing RF power on the other, the remaining RF power (linear damping on electrons and ions, if any) is adjusted self-consistently. Results from a GENRAY run are shown in Fig. 1. The antenna is shifted down by 1 m, which allows a better access of the plasma center by rays. In this GENRAY run, almost all RF power is deposited to electrons. Other species in ray-tracing simulations are thermal D + , T + , thermalized He-4, and a small fraction of Be-9 and Ar-40 that make up Z eff = 1.7. In all simulations discussed in this paper we use the total RF power of fast waves equal to 10 MW, and the initial parallel refractive index is set to n || = 1-5. For the rays shown in Fig. 1 the profile of power absorption by electrons is quite broad from plasma center to the edge, with a peak of power being at  0.63 ( is defined through the sqrt of toroidal magnetic flux) and another narrow peak present at  0.05.
The results from CQL3D runs are shown in Fig. 2. The profile of RF power absorbed by alphas in Fig. 2

(a)
is based on the solution of FP equation obtained by CQL3D-FOW, while the profile for deposition on electrons is based on self-consistent linear damping of the power flowing in rays, using the damping rate from GENRAY. It is seen that as the power deposition to alphas grows in time (thin lines converging to the bold line) the power absorbed by electrons is diminishedmainly at  < 0.6. The damping on alphas is mostly caused by resonance with their 16 th IC harmonic. The parasitic absorption by alphas is very high, 36% of the total injected RF power. This is related to penetration of waves into the plasma core where the content of alphas is relatively large. In a modified launching condition, when RF antenna is not shifted from the midplane, the rays tend to miss the core. As a result, in the GENRAY run there is no deposition of power to electrons at  < 0.2. From the CQL3D calculations, the damping on alphas is reduced from 36% to 30%.
A better result in power partition is achieved with higher wave frequency f RF = 800 MHz. The profiles of power deposited to alphas and electrons are shown in Fig. 3. The antenna is positioned at the midplane in this run. The power absorbed by alphas is 14%. If the antenna is shifted by 1 m down off the midplane, it becomes somewhat larger, 19% of the total RF power. The improvement is due to two factors. First, at higher frequency, the helicon waves tend to propagate more poloidally than radially, which moves the region of power deposition away from the plasma core. Second, and more important factor is that the waves interact with higher harmonic numbers: / cα = 21-28 along rays for the f RF = 800 MHz runs, which reduces the rate of alpha IC heating.
Another concern in parasitic power absorption is the damping on fast D + ions produced by NBI. In a separate run, using 33 MW NBI source, we obtained a 5% damping on D + . It is located at   0.65-0.95, and is caused by resonances at 24 th -27 th IC harmonics of D + .

ITER scenario EC_60s_1f
Different from the previous case, this scenario [5] has a lower auxiliary power of 56.5 MW (16.5 NB + 20 EC + 20 LH). Although the central magnetic field is unchanged, 5.3 T, the plasma current is lower, 9 MA. The central temperature is almost the same, T e0  T i0  36 keV, but the density is reduced to n e0 = 6.0e19 m -3 , n D = n T  2.4e19 m -3 . Because of the lower density of D + and T + , the power of fusion alphas is only P -source = 57.9 MW. One would expect a significantly smaller absorption of RF power on alphas, but it is not the case. The results from GENRAY run for a similar case of Fig. 1 with f RF = 500 MHz and the vertically shifted antenna are shown in Fig. 4. The corresponding profiles of power deposition obtained with CQL3D-FOW are shown in Fig. 5. The absorption of RF power on alphas is the same as in the previous scenario, i.e., 36%. This can be explained by observation that because of the lower density, the rays' radial group velocity is higher, and the power is deposited into the interior region with the larger population of alphas. The change of the power deposition profile is noticeable from comparison of Fig. 2 and Fig. 5. In the present case the peak of the p e (, t=0) profile is at  = 0.48 (Fig. 5(b), the top curve), while in the previous scenario the peak is at  = 0.63 ( Fig. 2(b)). The peak position of the power deposition to alphas is also shifted from  = 0.57 to  = 0.48 in the present scenario, and there is an additional absorption at  < 0.3, as seen in Fig. 5(a). For the alpha power absorption, the most of damping occurs at / cα = 16 resonance layer, which is shifted by 15 cm inward (  -0.1) for present equilibrium conditions, comparing to the previous scenario.
Switching to the higher frequency of f RF = 800 MHz yields a better power partition in favor of electrons. The results for the case of the midplane-centered antenna are shown in Fig. 6. The absorption on alphas is only 11%. This is the lowest fraction of RF power absorbed by alphas among the two ITER scenarios with different frequencies (500 or 800 MHz) and antenna positions. If the antenna is shifted by 1 m down off the midplane (in the present scenario with 800 MHz) the power deposited to alphas becomes 1.7 MW (17%). An additional CQL3D-FOW calculation for the power deposition to fast D + ions produced by 16.5 MW NBI source (half of that in the previous scenario), gives only a 1% of RF power damping on D + with two narrow peaks in the plasma edge region around   0.8 and   0.9. Thus, the result shown in Fig. 6 appears to be optimal for electron CD by the helicon waves. The profile of current based on GENRAY calculations is shown in Fig. 7. The CD efficiency is 0.047 A/W (total RF current is 0.47 MA), based on asymptotic formula by Ehst-Karney [7]. The linear CD calculation is sufficient for present purposes, since the electron quasilinear diffusion is weak. This CD efficiency is comparable to that expected for LHCD in ITER (0.04-0.06 A/W).
It should be mentioned though, that the distribution of alphas is generating a current as well, as they slow down and are transported in the radial direction. The profile of this current, derived from the solution of the FPE by CQL3D-FOW, is shown in Fig. 8. The value of this current reaches 0.37 MA at t = 4 sec (0.24 MA if the electron screening current is taken into account). It has nothing to do with the RF power, but it is rather a result    Returning to the optimal CD scenario with f RF = 800 MHz, 1.1 MW deposited to alphas, it is instructive to look at the distribution function of alphas. For the peak position of power deposition at the midplane coordinate R 0 = 779 cm ( = 0.65) the distribution function is shown in Fig. 9. It is the steady state solution of FPE at t = 2 sec (except the thermalized population of alphas at u/c < 5e-3 that keeps increasing). It is seen that the waves produce a high-energy tail beyond the birth energy of 3.52 MeV, especially at pitch angles close to π/2. At energies below 3.52 MeV, the distribution has a strong variation in pitch angle. This variation originates from the similar variation of the source operator: the copassing and especially co-going barely trapped particles have excursions into the hot inner region of plasma which is well populated by alphas. In contrast, the counter-going particles come from a relatively cold edge region. Consequently, the source operator is overpopulated at pitch angles less than π/2 (co-going side) and under-populated at pitch angles greater than π/2 [8]. This is another important feature (besides the neoclassical transport) that appears in FOW formulation.

Conclusion
Using high-harmonic or helicon fast waves for electron current drive in ITER is a promising CD method, but a big concern is the parasitic absorption by fusion-born alphas. The modeling with the CQL3D-FOW bounceaverage Fokker-Planck code for the selected two ITER scenarios proves that the absorption of RF power on alphas can be as high as 36%. To reduce this absorption, the RF frequency has to be increased to at least 800 MHz. In the best case that we could achieve the absorption on alphas was 11% of the total RF power. The vertical displacement of RF antenna can serve as an additional tool for changing in some extend the power partition between electrons and alphas, at a cost of changing the localization of the absorption power profile.
Each CQL3D-FOW simulation described in the paper is done in a relatively short computational time: 1.5 hours using 192 cores. The steady-state solutions for a physical time t = 2 sec are obtained in 20 time steps. Computer runs with the time step dt = 0.1 sec (and larger) are stable because of the fully-implicit solver. The code can be applied for other ITER scenarios in future work. In general, it is demonstrated that the code is a valuable tool for modeling of fast ions' RF and transport physics in ITER and other tokamaks.