Accelerated procedure to solve kinetic equation for neutral atoms in a hot plasma

The recombination of plasma charged components, electrons and ions of hydrogen isotopes, on the wall of a fusion reactor is a source of neutral molecules and atoms, recycling back into the plasma volume. Here neutral species participate, in particular, in charge-exchange (c-x) collisions with the plasma ions and, as a result, atoms of high energies with chaotically directed velocities are generated. Some fraction of these hot atoms hit the wall. Statistical Monte Carlo methods normally used to model c-x atoms are too time consuming for reasonably small level of accident errors and extensive parameter studies are problematic. By applying pass method to evaluate integrals from functions, including the ion velocity distribution, an iteration approach to solve one-dimensional kinetic equation [1], being alternative to Monte Carlo procedure, has been tremendously accelerated, at least by a factor of 30-50 [2]. Here this approach is developed further to solve the 2-D kinetic equation, applied to model the transport of c-x atoms in the vicinity of an opening in the wall, e.g., the entrance of the duct guiding to a diagnostic installation. This is necessary to determine firmly the energy spectrum of c-x atoms penetrating into the duct and to assess the erosion of the installation there. The results of kinetic modeling are compared with those obtained with the diffusion description for c-x atoms, being strictly relevant under plasma conditions of low temperature and high density, where the mean free path length between c-x collisions is much smaller than that till the atom ionization by electrons. It is demonstrated that the previous calculations [3], done with the diffusion approximation for c-x atoms, overestimate the erosion rate of Mo mirrors in a reactor by a factor of 3 compared to the result of the present kinetic study.


Introduction
In devices for thermonuclear fusion of hydrogen isotopes deuterium and tritium these are in the form of a hot strongly ionized plasma. To avoid too strong interactions of this plasma with the machine wall there is normally a peripheral edge layer, the so called scrape-off layer (SOL), where plasma particles and energy is transported along magnetic field lines into a special volume, the so called divertor [4]. In a future reactor like DEMO [5,6] processes in the divertor will significantly restrain this flow in the main part of the SOL [2] and a significant fraction of plasma particles, diffusing from the plasma core, will reach the vessel wall before they are exhausted into divertor. Arising plasma fluxes to the wall, with a density of 10 20 m −2 s [2], saturate the wall with fuel particles in a time much shorter than the discharge duration. Thus, a comparable influx of neutral species will recycle back into the plasma. These neutral particles are not confined by the magnetic field and can penetrate deep, at several centimeters into the SOL plasma. Here charge-exchange interactions with plasma ions generate atoms of energies much higher than that  Statistical Monte Carlo methods [7] are normally used to model c-x atoms at the edge of fusion devices. A crucial obstacle to apply these approaches routinely for extensive parameter studies, e.g., with the aim to optimize the duct geometry, is too large calculation time for a reasonably small level of accident errors. In a one-dimensional geometry an alternative approach based on iteration procedure to solve the kinetic equation represented in an integral form has been elaborated in Ref. [1]. Being, unlikely to Monte Carlo, free from statistical noise and permitting the convergence of iterations to the error level defined by the machine accuracy, this method, nonetheless, is also time consuming. However, by applying the pass method to evaluate integrals from functions, including the velocity distribution of plasma ions, a tremendous acceleration of kinetic calculations, at least by a factor of 30-50, has been achieved in Ref. [2]. This has allowed to perform a thorough comparison of calculations for 1-D transfer of c-x atoms in a hot plasma in the DEMO SOL done both with the kinetic description and by applying the so called diffusion approximation. Strictly speaking, the latter is relevant under plasma conditions of low temperature and high density, e.g., in the divertor, where the mean free path length between c-x collisions is much smaller than that till the atom ionization by electrons. Nonetheless, as it has been demonstrated in Ref. [2], even for the conditions of a hot SOL plasma the results of calculations done with the diffusion approximation deviate less than by 30% from those obtained kinetically and can be successfully applied for first estimates.
For diverse purposes numerous openings will be made in the reactor wall, e.g., for ducts leading to optical diagnostic installations collecting the radiation from the plasma. Such installations will be protected from the direct impact of charged plasma particles, moving mostly along the magnetic lines and penetrating into ducts on a Larmor radius depth of 1mm. Hot c-x atoms, unconfined by the magnetic field, can, however, freely hit and erode these installations. Since at the position of an opening the inflow of recycling neutrals, and, thus, the source of c-x atoms are absent, the transfer of the latter has two-or even three-dimensional pattern. A 2-D situation in the vicinity of a circular opening has been modeled in Ref. [3], by solving a 2-D diffusion equation for c-x atoms. In spite of the relevance of the diffusion approximation in the 1-D geometry, considered in Ref. [2], its applicability for multi-dimensional situations is, however, still questionable. In the present paper the approach developed in Ref. [2] to accelerate 1-D kinetic calculations is elaborated further to solve a 2-D kinetic equation for modeling of c-x atoms in the vicinity of a circular opening in the wall. The results of calculations are compared with those found in Ref. [3] by applying the diffusion approximation for c-x atoms. The remainder of the paper is organized as follows. In the next section the basic equations are formulated and reduced to the form used for the numerical treatment. A pass method applied to asses the integral in the velocity space is also discussed here. In section 3 the results of calculations both with kinetic and diffusion description of c-x atoms are presented and compared. Conclusions are drawn in section 4.

Basic equations
Henceforth we operate in the coordinate system shown in figure 1 and characterized by the distances x and ρ from the wall and opening axis, respectively, under the assumption of the opening axis-symmetry. The velocity distribution function of c-x atoms, f (x, ρ, V x , V ρ ), is governed by the kinetic equation: where ν a = n (k a ion + k a cx )is the total frequency of the atom disintegration in the plasma through the processes of ionization by electrons and charge-exchange with ions, characterized by the rate coefficients k a ion and k a cx , respectively; n is the plasma density, the same for electrons and ions due to quasi-neutrality; S cx = S 0 cx + nk a cx n cx is the density for c-x atoms source where the former contribution S 0 cx , specified in Refs. [2,3], is due to charge-exchange of recycling neutrals, and the latter is owing to further c-x interactions of atoms generated from recycling species; the plasma ions with the mass m i and temperature T i have a Maxwellian velocity distribution with the thermal velocity V th = 2T i /m i .

Reduction to integral-differential equations
Henceforth, V ρ is discretized as ±V l ρ , with V l ρ = ΔV (l − 1/2), ΔV = V max /l max and l = 1, · · · , l max . The distribution functions ϕ ± l of particles with ρ-velocity component in the range are governed by the equations following from the integration of equation (1) over these ranges: where S cx,l = S 0 cx + nk a cx n cx δ l , with being the fraction of ions in the the l-range.
From equations (2) one gets equations for the new variables The latter equation is straightforwardly integrated with respect to x, providing x 0 being the position where the boundary conditions are imposed: since there is no influx of c-x atoms from the wall, γ l (x 0 = 0) = 0 for V x > 0; besides neutral species are absent deep enough in the plasma, i.e. γ l (x 0 → ∞) → 0 for V x < 0. By substituting the found expression for γ l into equation (3), one gets: In the second term on the right hand of the equation above both factor are expanded into the Taylor series up to linear terms, . By taking into account that that the width of the region occupied by neutral atoms is of several mean free path lengths (MFPL), i.e. typically |x − . As a result we get the following equation: This can be integrated as the first order equation with respect to x, with the same boundary conditions as for γ l , i.e. ϕ l (x = 0) = 0 for V x > 0 and ϕ l (x → ∞) → 0. Finally, for the density of atoms with one obtains the following integral-differential equation: where The total density of c-x atoms n cx = lm l=1 η l . Equation (7) is reducing to an integral one obtained in the 1-D case, either equation (16) in Ref. [1] or equation (11) in Ref. [2], by neglecting the ρ-derivative and considering one V l ρ -value only, i.e. by assuming l max = 1. In this case l ≡ 1, δ l = 1, n cx = η 1 and S cx,1 = S 0 cx + nk a cx n cx . As it was demonstrated in Ref. [2] this integral equation can be reduced to a ordinary differential diffusion equation if (i) the electron temperature is low and the ionization rate is much smaller with D th = V 2 th / (2ν a ). The boundary conditions to equation (8) imply: (i) c-x atoms leave the plasma with the velocity close to the ion thermal one; (ii) neutral species are absent far from the wall, n cx (x → ∞) = 0. Additionally, ∂ ρ η l , ∂ ρ n cx = 0 for ρ = 0 and ρ → ∞ for equations (7) and (8) .

Assessment of the velocity space integral I α
By inspecting equation (7) one can see the origin of large calculation time for kinetic computations of c-x atoms. Even in the 1-D case for each grid point one has to calculate enclosed double integrals over the ion velocity space and over the whole plasma volume, 0 ≤ y ≤ r w , and repeat this in an iterative procedure with respect to the whole profiles of n cx . For the 2-D situation the runs through the grid points in the (ρ, V ρ )phase subspace are involved into calculations additionally; however these can be efficiently parallelized. Fortunately, the integral I α , resulting from the generation of c-x atoms with the ion Maxwellian distribution function and the consequent attenuation in the ionization and c-x collisions, can be assessed with a high enough accuracy by an approximate pass method [8]. The function F α (u) is shown in figure  2a for α = 0.3, 1 and 3. One can distinguish four u-intervals where F α (u) behaves principally differently: u ≤ u 1 , u 1 < u ≤ u m , u m < u ≤ u 2 and u 2 < u, with the points u 1,2 and u m defined by the conditions F α (u m ) = 0 and F α (u 1,2 ) = 0 for the derivatives F α = −F α g 1 /u 2 , F α = F α g 2 /u 4 where g 1 (u) = 2u 3 + u − α, g 2 (u) = (g 1 + u) 2 − 2 + 6u 2 u 2 , u m is a root of a cubic equation, u m = 3 √ b + a − 3 √ b − a with a = α/4 and b = a 2 + 1/216. By searching for u 1,2 , one can see that for the roots of interest the equation g 2 (u 1,2 ) = 0 reduces to the following ones: h 2 (u 2 ) = g 1 (u 2 ) + u 2 − u 2 2 + 6u 2 2 = 0. By selecting properly the initial approximations to u 1,2 , these equations are solved with the Newton method, providing accurate enough solutions after 3-5 iterations.
To assess the integrals I α we approximate function F α (u) as a linear one at u ≤ u 1 and by polynomials of the fifth order at the intervals u 1 < u ≤ u m and u m < u ≤ u 2 with the coefficients selected to reproduce F α,m,1, This results in: where δ 1,2 = |u m − u 1,2 |. Figure 2b shows I α versus α found with the pass method outlined above and evaluated numerically. By approximating the direct numerical results very accurately, the procedure above allows to reduce the CPU time by a factor of 30-50 at least. As an alternative to the usage of formula (9) one can calculate in advance the integrals I α (α) and save this in a table. Then, for any particular α, the integral I α (α) can be interpolated from the values in the table for arguments larger and smaller than that in question. However, the searching of the correct α-interval in this table requires also some time which grows up with α → 0 since I α → ∞ as E 1 (α) and the density of data in the table has to be increasing. Thus, it is not obvious that the usage of tables is more time saving than formula (9).

Results of calculations
The radial profiles of the plasma parameters used for calculations of c-x atom parameters are presented in figure 3 and were computed in Ref. [2] with the input parameters from the European DEMO project [5,6].
In figure 4 we demonstrate the variation of the c-x atom density with the distance from the wall, far from a circular opening of the radius ρ 0 = 0.05m and at the opening axis. There is a noticeable difference between the profiles calculated kinetically and in the diffusion approximation.
To comprehend the importance of this difference we assess the rate of erosion by c-x atoms of a Mo mirror imposed in the duct at the distance h from its opening. In a strict diffusion approximation the energy distribution function of c-x atoms coincides in any point in the plasma with that of the ions. Thus, the averaged energy of atoms entering the opening in the wall is close to the ion temperature at this position, T i (0). However, the ion temperature increases by several times at the penetration depth of c-x atoms, compare figures 3 and 4. C-x atoms generated here partly escape without further collisions with electrons and ions. Therefore, the distribution over the energies E of atoms, hitting the opening, can significantly differ from the local Maxwellian one of magnetized ions: Consider an infinitesimally thin ring of the width dρ, thickness dx, with the radius ρ and situated in the plasma at the distance x from the wall. C-x atoms of the energy within the range E, E+dE are generated in this ring at the rate: Since c-x atoms move in all directions some of them hit the mirror inside the duct, see figure 1.
Henceforth we consider the surface position at the duct axis only. A simple geometric analysis shows that c-x atoms from the rings with ρ ≤ ρ max = ρ 0 (1 + x/h) can hit this surface point. The contribution of the plasma volume in question to the density of the c-x atom flux perpendicular to the installation surface is: where s = 1 + ρ 2 / (x + h) 2 −1/2 is the cosine of the atom incidence angle with respect to the port axis. The exponential factor in equation (12), with λ = x 0 ν a dx / 2E/m, takes into account the destruction, due to ionization and charge-exchange collisions, of the atoms in question on their way through the plasma. By ionization the atom completely disappears. By charge-exchange a new atom is generated. The latter process is taken into account by the second term in the source density, S cx = S 0 cx + nk a cx n cx , and by the summation/integration of dj cx contributions from rings situated in the plasma at all possible x.
The energy spectrum of atoms hitting the mirror is characterized by their flux density in the energy range dE, γ = dj cx /dE, where the integration is performed over the ring radius ρ and distance from the wall x. By proceeding from ρ to s, according to the relation ρ = (x + h) 1/s 2 − 1, one gets: The density of the outflow of the mirror material particles eroded by physical sputtering with c-x fuel atoms can be assessed as follows: Here Y sp is the sputtering yield, whose dependence on the projectile energy E and the cosine s of the incidence angle is calculated by applying semi-empirical formulas from Ref. [10]. The erosion rate h sp = Γ sp /n sp , with n sp being the particle density of a mirror of Mo and measured in nm/year, is demonstrated in figure 5 as a function of the installation surface position h inside the port, for n cx (x, ρ) involved into S cx calculated both kinetically from equations (7) and in the diffusion approximation, by integrating equation (8). rate found with kinetically assessed n cx (x, ρ) is noticeably below that determined by using the diffusive approximation for c-x atoms in the plasma. The iteration approach elaborated in Ref. [1] to solve 1-D kinetic equation for c-x atoms has been developed further to describe the transport of these species in 2-D situation, in particular, in the vicinity of a circular opening in the wall of a fusion reactor. Unlike to Monte Carlo methods the iteration one does not generate statistical errors. To perform a thorough comparison with the results obtained by using a diffusion approximation for c-x atoms [3], the solving procedure is accelerated by a factor of 50, by applying pass methods to assess integrals in the velocity space [2]. Calculations done by both approaches demonstrate that the diffusion approximation provides significantly overestimated estimates for the erosion rate of a first mirror in a fusion reactor like DEMO compared to the kinetic calculations. Thus, in spite of the relevance of the diffusion approach for c-x atoms in 1-D calculations of the plasma parameters in the SOL, see Ref. [2], the kinetic one seems to be unavoidable for calculations of characteristics being very sensitive to the c-x atom energy like the erosion rate of the reactor wall and diagnostic installations.