A full wave solver integrated with a Fokker–Planck code for optimizing ion heating with ICRF waves for the ITER deuterium–tritium plasma

Efficient ion heating is crucial for future fusion devices, and the only way to heat ions directly is ion cyclotron resonance heating. Reported here is a full wave solver integrated with a Fokker–Planck code for optimizing ion heating with ion cyclotron range of frequency waves for the International Thermonuclear Experimental Reactor deuterium–tritium plasma. Both the direct absorption of minority ions and the power transfer to bulk ions via collisions are considered, while also accounting for the edge effects on ion absorption near the core. The simulation results show that the appropriate scrape-off layer density profile and parallel wave number lead to enhanced edge coupling and broaden the absorption region with moderate absorption intensity of the minority ions, which is very important for ion heating. More power from the heated ions is transferred to bulk ions than to electrons through collisions in our simulation via optimization, and reducing the total RF power results in a significant increase of the absorbed fraction of bulk ions.

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Efficient ion heating is crucial for future fusion devices because it offers direct increase of ion temperature and increased fusion yield [1][2][3], and high-Q operation can be obtained faster than with electron heating.Ion cyclotron resonance heating (ICRH) is the only method that can heat ions directly via the ion cyclotron resonance mechanism between ions and waves.Moreover, fast waves in the ion cyclotron range of frequency (ICRF) can propagate into the core of a highdensity plasma without a high-density cut-off [4], thereby facilitating effective ion heating.ICRH is one of the three main auxiliary heating methods in the International Thermonuclear Experimental Reactor (ITER) [5], and the planned ICRF antenna in ITER is designed to provide a radiated power of 5-20 MW in the frequency range of 40-55 MHz [6][7][8][9][10].Secondharmonic heating of tritium ions is considered to be the main ICRF heating scheme in ITER in the burning plasma, but during the initial ramp-up phase, the plasma temperature and density are rather low, resulting in weak harmonic heating, and so the minority-ion heating scenario becomes dominant [11][12][13][14][15].A small amount of 3 He ions is puffed into the deuteriumtritium (D-T) plasma, and these minority ions absorb the radio frequency (RF) power via the fundamental resonance and then transfer energy to the D-T bulk ions via Coulomb collisions.However, the RF power can also be absorbed directly by electrons via a mechanism such as Landau damping.Besides, the energy of the heated minority ions is also transferred to electrons through slowing down with Coulomb collisions.In ICRH, there is a competition between ions and electrons regarding wave power absorption.Therefore, optimizing ion heating in the minority-ion heating scenario involves two aspects: (i) increasing the absorption of wave power by minority ions rather than electrons and (ii) increasing the power fraction that the heated minority ions transfer to the bulk ions via collisions.If the average energy of the heated ions is lower than the critical energy then more power is transferred to the ions, which implies that moderate locally absorbed power density and broad absorption region of the minority ions should benefit bulk-ion heating.
In the initial ramp-up phase for ITER with relatively low plasma density and temperature, previous simulation results have shown that the bulk-ion absorption fraction is ca.50% in the minority-ion heating scenario [16].A novel ICRF bulkion heating method in D-T plasma using impurities such as Be or Li has been proposed [17,18]; the critical energy of Be ions is almost three times that of He ions because of the higher atomic mass of Be, and this new three-ion heating scenario is believed to have the potential to provide a larger fraction of fuel ion heating.However, it has recently been planned to replace Be with W as the first wall material in ITER, so it is meaningful to carry out further research on optimizing the ion heating in the D-T-( 3 He) minority-ion heating scenario.In addition, the numerical simulations of ICRF heating for ITER are commonly performed by core solvers.The twodimensional (2D) full-wave codes TORIC [19,20] and EVE [21] have been used to simulate the D-T-( 9 Be) minority-ion heating scenario in ITER, accurately calculating the absorption coefficients and power redistribution of the heated ions.Various ICRH scenarios for ITER have been evaluated by PRETOR + PION [22][23][24][25][26], and the absorbed power of various species for the D-T- 3 He-Be-α- 4 He scenario has been shown in detail [3].However, the scrape-off layer (SOL) plasma model and RF antenna structures in these core solvers are too simple, and the edge coupling cannot be analysed well by these core codes.Instead, the effects of poloidal phasing on ICRH power absorption and coupling for ITER are studied using a novel code (FEMIC) based on the finite-element method that includes both core and edge models [27].We applied the finiteelement method based on the approach of [27] to simulate ICRH on EAST [28,29].
Herein, we report the development of the aforementioned wave solver integrated with a quasi-linear Fokker-Planck code for simulating the optimization of ICRF ion heating in the ITER D-T plasma including the effect of edge coupling on core ion heating, and both direct absorption of minority ions and the power transfer to bulk ions via collisions are considered.This paper is arranged as follows: in section 2, the physical model is described; in section 3, the simulation results are presented and the minority-ion heating is discussed; in section 4, results are given for power redistribution of heated ions; finally, some conclusions are drawn in section 5.

Full-wave model
The quasi-uniform approximation is used to simplify the solution of the full-wave equation, as described in [27], and this approach involves pre-calculating the wave vectors to obtain the local dielectric tensor.The assumption of toroidal symmetry is applied.We Fourier decompose in the toroidal direction and calculate the wave field for a single toroidal mode number [21], and the propagation and absorption of waves are simulated on the 2D poloidal cross-section of the ITER tokamak.The full wave equation can be expressed as where E is the electric field of the wave, ω is the wave's frequency, J ant is the current of the antenna, and k ⊥ , k // are the wave vectors in the vertical and parallel directions, respectively.The term K is the dielectric tensor based on the Maxwellian distribution and given by where I is the identity matrix and χ j is the hot susceptibility with Maxwellian distribution for species j; the detailed expressions for hot susceptibility are given in [12].The FEMIC allows the wave equation to be solved, and the wave propagation and absorption in both the SOL and core are described well.

Quasi-linear Fokker-Planck model
Following Brambilla [30], we solve the quasi-linear kinetic equation on each magnetic surface with cubic finite elements.
The toroidal trapping effect of fast ions on ICRH is neglected.The quasi-linear diffusion coefficient (surface-averaged) in the diffusion operator is constructed by inputting the wave electric field and the absorbed power in the resonance region, which are the outputs of the full-wave code, see formula (6) for details.Meanwhile, a linear collision operator is adopted in the equation, assuming that the velocity distribution of fast ions maintains a steady state by dissipating energy on the background ions and electrons.The Fokker-Planck equation can be written as 0 = where F i is the velocity distribution function of the test ions.
The first term on the right-hand side of the equation is the quasi-linear diffusion operator that describes the ion cyclotron resonant interaction between fast waves and test ions.The second term is the linear collision operator that describes the Coulomb collisions between the fast ions and the bulk plasma, and S i is a source term.The quasi-linear diffusion operator [30] is and the quasi-linear diffusion coefficient is expressed as where E + is the left-handed electric field of the wave, which plays an important role in ion cyclotron resonance, J is the Bessel function (for the p th harmonic frequency, where p = 0 indicates the fundamental resonance), m and n are respectively the toroidal and poloidal mode numbers of the Fourier expansion, where λ P = and ω = v ⊥ /v thi , D 0 can be evaluated according to energy conservation between the surface-averaged absorbed power from the wave solver and the heating rate [31].Then the quasilinear operator is represented in spherical coordinates (7) The collision operator is expressed as where v is the velocity of the test ions, µ is the cosine of the pitch angle, and the expressions of the collision coefficients as Ψ c (v), Ψ τ (v), Θ c (v) can be found in [32,33].The distribution function F i (v, µ) is represented as an expansion in Legendre polynomials regarding the velocity pitch angle, typically incorporating five to eight terms The steady-state solution is reduced to the solution of a set of coupled ordinary differential equations, which can be solved by the cubic Hermite finite elements method [34].Due to the quasi-linear kinetic equation is solved on each magnetic surface, the solved ion distribution function is in 2D velocity space and 1D radial space.

Magnetic equilibrium, geometry and plasma parameter
We focus on the ramp-up phase of ITER with relatively low density and temperature compared to the burning phase, while the rate of D-T fusion reactions is small.The adopted EFIT equilibrium document, plasma temperature and density in our simulation are from the IMAS database [35].#120012 L-mode for ITER.The toroidal magnetic field is B 0 ∼5.3 T, the central density is n e0 = 4.61 × 10 19 m −3 , the central electron temperature is T e0 = 5.3 keV, and the central ion temperature is T e0 = 5.0 keV.The ICRH power is normalized to 1 MW in the full wave solver.The geometry of the poloidal cross-section of ITER is shown in figure 1, and the 2D profiles of plasma parameters are shown in figure 2. The first wall, divertor and antenna box are represented by black lines and are assumed to have perfect electrical conductivity, and the magnetic surfaces are represented by colour lines.The antenna straps are represented by green lines and are assumed to be infinitely thin metal plates (consisting of six poloidal radiating straps) [35].The current density in current straps is given as Where t ant is the unit tangential vector, Jk is the current density amplitude and θ k is the poloidal phasing of the antenna.For simplicity, only the main radiation conductor is retained in the simulation model, and the Faraday screen is not considered.Assuming that the current straps are placed in a vacuum environment, the grid is configured as a triangular grid with a maximum size of less than 1 cm.The impact of edge coupling on core ion absorption can be taken into account in this model.The model requires ca.80 GB of RAM, and the average simulation time for a single run is ca.17 min.

Benchmark
Vallejos used the FEMIC to study how phase affects ion cyclotron heating in the ITER device [27], while Zhang further verified its effectiveness by simulating ion cyclotron heating in the EAST device [28,29].To ensure the reliability of the simulation results in this paper, benchmarking is carried out between the FEMIC code integrated with the Fokker-Planck code (referred to herein as FEMIC-FP) and the TORIC-SSFPQL code for ICRH in ITER.The SOL was neglected in both codes in benchmark, antenna frequency is 53 MHz, ICRH power is normalized to 1 MW, the toroidal mode number is 33, the 3 He ion concentration is 5%, and the D:T ratio is 1:1.
The power absorption fractions of various species with different wave frequencies simulated by FEMIC and TORIC are compared in figure 3, the power density profiles calculated by the two codes are shown in figures 4(a) and (b), the collisional energy redistribution and profiles of collisional energy transfer calculated by FEMIC-FP and TORIC-SSFPQL can be seen in table 1 and figure 6, and the comparison of the 2D velocity distribution of minority ions is shown in figure 5.The results simulated by FEMIC-FP and TORIC-SSFPQL agree well, and the discrepancies can be explained by the difference in the plasma models.

Effect of wave frequency on minority ion absorption
For a given magnetic configuration, the choice of wave frequency determines the location of the ion cyclotron resonance layer.Figure 7 shows the absorption faction with wave frequency ranging from 38 MHz to 65 MHz.In the wave frequency range of 45-56 MHz, the absorption fraction of minority ions is close to 80%.The power deposition in the poloidal cross-section for a frequency of 53 MHz is shown in figure 8(a), which indicates strong minority-ion absorption with a broadened resonance layer near the core; the electron damping region is wide, but with weak absorption.For comparison, at a wave frequency of 43 MHz as shown in figure 8(b), the ion cyclotron resonance layer is located offaxis on the low-field side, and the absorption of wave power by ions and electrons is almost equal.Therefore, the wave frequency of 53 MHz with a central resonance is adopted in our simulations.

Effects of SOL density profile on coupling and absorption
Fast waves emitted from the ICRH antennas on the low-field side first encounter the right-hand cut-off layer in the SOL, and tunnel though the evanescent region into the plasma, so SOL density profile obviously affects the edge coupling and also influences the core ion absorption of the wave.The FEM code includes the SOL and core plasma and so can be used to study how the SOL density profile affects the edge coupling and core ion absorption.We use the coupling power P c to describe the coupling capability of ICRF waves, where P c is the total coupled power per unit current density.The electron density profile in the simulation is expressed analytically as where n e is the electron density in the core, n eL is the electron density at the last closed flux surface, n min is the minimal electron density within the SOL, and L is the decay length of density within the SOL.First, we simulate the effects of minimum density on edge coupling.As can be seen in figure 9, the coupling power P c increases gradually with increasing minimum density, which means that the coupling effect is better at a larger minimum density in the SOL.This is consistent with the good coupling effect obtained in experiments such as JET [36] and AUG [37,38] by increasing the edge density by injecting gas near the antenna.However, increasing the minimum density at the edge can lead to possible damage to the antenna and the first wall.Therefore, the effects of the SOL density gradient on coupling and core absorption are investigated next.We can modify the density profile of the SOL plasma by changing the decay length L (see equation ( 11)), where L is inversely proportional to the density gradient, as shown in figure 10(a).The effect of the density gradient in the SOL on the coupling power P c is shown in figure 10(b).As L increases, the density gradient decreases and P c increases gradually, which means that the coupling effect is enhanced with a small density gradient.The evanescent length becomes short with a small density gradient (see figure 10(a)), which means that the wave can tunnel through the evanescent region into the plasma more easily.On the other hand, the wave may be obstructed and reflected when it encounters a large density gradient, resulting in a poor coupling effect.Figure 11(a) shows the left-handed electric field of the wave in the poloidal cross-section with a large density gradient (L = 0.01), the fast wave partially transforms into a kind of surface wave propagating along the outermost closed magnetic surface when it encounters the large density gradient, while the electric field of the wave propagating into the

Effects of parallel wave number on edge coupling and core absorption
The spectrum of the parallel wave number k // excited by the antenna is a key factor influencing wave propagation and absorption.The coupling power P c and the absorption fraction of the minority ions with different k // are shown in figure 13.P c decreases as k // increases because k // is proportional to the length of the evanescent layer, and the evanescent length is shorter with lower k // , and so the wave coupling effect is better.On the other hand, as k // increases, the absorption fraction of minority ions increases initially and then stabilizes at ca. 90% after k // reaches 5.This is because the Doppler broadening of the ion absorption region increases as k // increases, and so the ion absorption fraction increases.The peak of k // with several typical phases of the ITER antenna is shown in figure 13, and k // = 4.8 with antenna phase (0, π, π, 0) might be the optimal choice for ion heating considering both wave coupling and ion absorption.Figure 14 shows the power deposition on minority ions in the poloidal cross-section with the main k // of three typical phases in ITER, and the ion absorption region with k // = 4.8 is the largest among the three cases.

Effect of minority-ion concentration on ion absorption
In the minority-ion heating scenario, the concentration of minority ions significantly influences their absorption.Figure 15(a) shows how the power absorption fraction of various species depends on the minority-ion concentration.The absorption fraction of 3 He ions decreases in general as the minority-ion concentration increases, and the absorption fraction of minority ions exceeds 90% when the minority-ion concentration is 1%-2%.As shown in figure 15(b), the ion-ion hybrid resonance (IIHR) layer moves towards the high-field side as the 3 He ion concentration increases, the distance between the IIHR layer and the ion cyclotron resonant layer increases gradually, and the left-hand electric field of the wave weakens in the resonance region, resulting in a reduction in the fraction of ion absorption.

Optimizing bulk ion heating
Direct heating of minority ions was studied in the previous section, and the power redistribution of heated ions is discussed in this section.For optimal ion heating, the wave frequency is set as 53 MHz, the minority-ion concentration is 2%, the decay length in the SOL density profile is set as 0.05 for a small density gradient, and k // = 4.8 is used in the following simulation.The left-hand electric field of the wave and absorption power density of 3 He ions calculated by the full-wave solution are used as inputs for the Fokker-Planck code.The coupling calculation of the full wave solver and Fokker-Planck solver is done only once in this work, and we are studying that the dielectric tensor of the full-wave solution to incorporate non-Maxwellian distribution, to ensure iterative self-consistency between the full-wave code and the Fokker-Planck code in the future.For 20-MW RF power, the heated minority-ion velocity distributions with different pitch angles are shown in figure 16, and the effective perpendicular tail temperature is 123.3 keV, which is below the critical energy.As shown in table 2, ca.70% of the power is transferred to the bulk ions and 25% is transferred to electrons via collisions.About 67% eventually ends up in bulk D-T ions and nearly 29% in electrons considering direct power absorption.
Moreover, decreasing the total RF power results in a significant increase in the power fraction transferred to the bulk ions.As shown in table 3, for the case of P RF = 5 MW, ca.90% of the power is transferred to the bulk ions, and ca.85% eventually ends up in the bulk D-T ions.Reducing the RF power may have advantageous effects on bulk-ion heating during the Figure 14.Power absorption of 3 He ions in poloidal cross-section.The power fraction of direct absorption of 3 He ions is merely 82.4% for 10% minority-ion concentration, and the power fraction of 3 He minority-ion transfer to bulk ions is 70.7% in this case, which is lower than the power fraction transfer to bulk ions with 2% minority-ion concentration.Minority ions transfer more energy to themselves for the high-concentration case.

Conclusions
As reported here, the FEM code integrated with a quasi-linear Fokker-Planck code was developed to study the optimization of ion heating with ICRF waves for ITER.We focused on the D-T-( 3 He) minority-ion heating scenario in the initial rampup phase.First, for optimizing the minority-ion heating, the wave frequency was set as 53 MHz, and the effects of the SOL density profile on the coupling and absorption were discussed in detail.Smaller SOL density gradient leads to enhanced edge coupling and a broadened absorption region with moderate absorption intensity, which is very important for ion heating.On the other hand, a surface wave is generated when ICRF waves encounter a large SOL density gradient, which may induce potential damage on the wall.The effects of k // on wave coupling and ion absorption were both taken into account, and the antenna phase (0ππ0) with k // = 4.8 was considered as the optimal choice.The absorption fraction of minority ions exceeded 90% when the minority-ion concentration was 1%-2%.
Next, the full-wave solution was used as the input for the Fokker-Planck code, and the power redistribution of heated ions was discussed.For 20-MW RF power, the simulation results showed that more power of heated ions is transferred to bulk ions than electrons, and decreasing the coupled RF power resulted in a significant increase of absorbed fraction of bulk ions, with the fraction of ion absorption approaching nearly 90% for P RF = 5 MW.Finally, the power redistribution with different minority-ion concentrations of 2% and 10% were compared.In our study on optimizing ion heating, both the direct absorption of minority ions and the power transfer to bulk ions via collisions were taken into consideration, while also accounting for edge effects on core absorption.Our research offers a comprehensive scheme and parameter reference for future ITER-ICRH experiments that can be extrapolated to future fusion reactors.Full wave tools which consider realistic 3D geometries were developed recently [39], This 2D wave solver has the potential to be extended to the 3D case, more toroidal effect will be considered in the future.

Figure 1 .
Figure 1.Geometry of the ITER cross section.

Figure 3 .
Figure 3. Power fractions of various species as a function of frequency.Solid lines: results calculated by FEMIC code; dashed lines: results calculated by TORIC code.

Figure 4 .
Figure 4. Power absorption profile calculated by (a) FEMIC code and (b) TORIC code.The ICRH power is normalized to 1 MW.

Figure 7 .
Figure 7. Power absorbed faction of various species with different frequency.

Figure 9 .
Figure 9. Coupling power Pc as a function of n min .

Figure 10 .
Figure 10.(a) Radial density profile of SOL for different L, where the red arrow points to the cut-off layer density, and the black dotted and solid arrows correspond to the decay length of L = 0.01 and 0.05, respectively.(b) Coupling power Pc as function of L.

L. Yin et alFigure 11 .
Figure 11.Left-handed electric field of wave (a) and power deposition of 3 He ions (b) in poloidal cross-section with L = 0.01.

L. Yin et alFigure 12 .
Figure 12.Left-handed electric field of wave (a) and power deposition of 3 He ions (b) in poloidal cross-section with L = 0.05.

Figure 13 .
Figure 13.Coupling power Pc, 3 He absorption (left axis) and 3 He absorption fraction (right axis) as functions of k // .The hollow circles indicate the 3 He absorption corresponding to the main k // of the typical toroidal phasing case.

L. Yin et alFigure 15 .
Figure 15.(a) Power absorption fraction of various species as 3 He concentration increases from 1% to 20%.(b) Positions of IIRH layer and ion cyclotron resonance layer.

Figure 17 . 3
Figure 17. 3 He ion velocity distribution for coupling RF power of 5 MW.

Table 2 .
Absorbed power fraction of ions and electrons for 3 (He)-D-T heating scheme with X 3 He = 2% and P ICRF = 20 MW.

Table 3 .
Absorbed power fraction of ions and electrons with different P ICRF .

Table 4 .
Absorbed power fraction of ions and electrons with 2% 3 He and 10% 3 He.