Noise-Free Rapid Approach to Solve Kinetic Equations for Hot Atoms in Fusion Plasmas Noise-Free Rapid Approach to Solve Kinetic Equations for Hot Atoms in Fusion Plasmas

At the first wall of a fusion reactor, charged plasma particles are recombined into neutral molecules and atoms recycling back into the plasma volume where charge exchange ( cx ) with ions. As a result hot atoms with chaotically directed velocities are generated which can strike and erode the wall. An approach to solve the kinetic equation in integral form for cx atoms, being alternative to statistical Monte Carlo methods, has been speeded up by a factor of 50, by applying an approximate pass method to evaluate integrals, involving the ion velocity distribution function. It is applied to two-dimensional transfer of cx atoms near the entrance of a duct, guiding to the first mirror for optical observations. The energy spectrum of hot cx atoms, escaping into the duct, is calculated and the mirror erosion rate is assessed. Computations are done for a molybdenum first mirror under plasma conditions expected in the fusion reactor DEMO. Kinetic modeling results are compared with those found with a diffusion approximation valid in very cold and dense plasmas. For ducts at the torus outboard a more rigorous kinetic consideration predicts an erosion rate by a factor up to 2 larger than the diffusion approximation.


Introduction
In devices for thermonuclear fusion research, for example, of the Tokamak type, particles of hydrogen isotopes, deuterium and tritium, are in the form of a hot fully ionized plasma [1,2]. To avoid a destruction of the machine wall, a special region, the so-called scrape-off layer (SOL), is arranged at the plasma edge, where particles stream along the magnetic field to the special target plates [3]. Normally, it is done by using additional magnetic coils to form a divertor configuration (see Figure 1a).
By reaching the divertor target plates, plasma electrons and ions are recombined into neutral atoms and molecules which are finally exhausted from the device by pumps. However, in a future fusion reactor like DEMO [1,2], only a minor fraction of 1% of neutrals generated at the targets will be pumped out. The rest of them is ionized again in the plasma near the targets. This "recycling" process significantly restrains the parallel plasma flow in the SOL [4]. Therefore, a considerable fraction of plasma particles lost from the plasma core will reach the vessel wall before they are exhausted into the divertor. Plasma fluxes to the wall saturate it with fuel particles in a time much shorter than the discharge duration, and a comparable amount of neutral species will recycle back from the wall into the plasma. Recycling neutrals are not confined by the magnetic field and penetrate at several centimeters into the SOL. Here, charge exchange (cx) collisions of them with ions generate atoms of energies much higher than that of primary recycling neutral particles. A noticeable fraction of such secondary cx atoms hit the vessel wall and erode it.
Statistical Monte Carlo methods [5] are normally used to model cx atoms at the edge of fusion devices. A crucial obstacle to apply these approaches for extensive parameter studies, for example, with the aim to optimize the duct geometry, has too long calculations needed to achieve reasonably small accident errors. This is, however, necessary, for example, to couple neutral parameters with the plasma calculations. 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 decades ago [6]. Being 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. The reason is the necessity to assess integrals in the velocity space from functions involving the ion velocity distribution function, additionally to integrations in the normal space. Recently [4], an approximate pass method has been applied to evaluate these integrals, and the acceleration of kinetic calculations by a factor of 50 has been achieved.
The amendments, outlined above, have allowed to perform calculations of the plasma parameters in the DEMO SOL with cx atoms described kinetically, by varying the input parameters, for example, the plasma transport characteristics, in a broad range [4]. In addition the results of these computations have been thoroughly compared with those obtained with cx species described in the so-called diffusion approximation. This approximation, often used in diverse edge modeling approaches to save CPU time (see, e.g., [7,8]), is strictly valid under plasma conditions of low temperature and high density where the time between cx collisions of atoms with ions is much smaller than that till their ionization by electrons. Across the DEMO SOL, the plasma density and the temperatures of electrons and ions change, however, by orders of magnitude [4].
The usage of a diffusion approximation for cx atoms, generated from species recycling from the wall, becomes especially questionable by considering the situation near an opening in the vessel wall. Such openings will be made in a reactor for diverse purposes, for example, for ducts leading to the first mirrors, collecting light emitted by impurity species in the plasma (see Figure 1b). These installations are inaccessible for charged plasma particles, moving mostly along magnetic field lines and penetrating into the duct at a distance of 1 mm. However, hot cx atoms, unconfined by the magnetic field, can freely hit and erode these installations. The erosion rate is very sensitive to the energy spectrum of cx atoms which can be significantly dependent on the modeling approach. At the position of the duct opening, the inflow of recycling neutrals is actually absent, and the transfer of cx atoms has two-or even threedimensional pattern. By calculating the density of cx atoms in the vicinity of a circular opening from a 2D diffusion equation, the erosion rate of a first mirror of Mo has been assessed in [9]. In the present paper, we extend the approaches, elaborated in [4] to model cx atoms kinetically in the 1D case, on a 2D geometry.
The results are compared with those of [9]. Although there cx atoms have been considered in the diffusion approximation to reduce CPU time, the new approach allows to perform more exact kinetic calculations even by orders of magnitude faster.

Basic equations
Although the concept of magnetic fusion is based on the idea that charged particles can be infinitely long confined within closed magnetic surfaces, there are diverse physical mechanisms leading to the losses across these surfaces (see, e.g., [10]). Therefore, electrons and ions escape through the separatrix into the SOL and may reach the first wall of the machine vessel.
Here, these species fast recombine into neutral atoms, with an energy comparable with that of ions. Partly, these atoms escape immediately back into the plasma, and these species are referred as backscattered (bs) atoms. The rest transmits the energy to the wall particles and may be bounded into molecules, desorbing back into the plasma as the wall, and is saturated with the working gas. Here, the densities of recycling backscattered atoms and desorbed molecules decay with the distance from the wall because diverse elementary processes such as ionization, charge exchange and dissociation in the latter case, caused by collisions with the plasma electrons and ions, lead to disintegration of neutrals [11]. By the dissociation and ionization of molecules, the so-called Franck-Condon fc atoms are generated of a characteristic energy E fc of 3:5 eV [3,11]. By charge exchange collisions of the primary neutral species, listed above, with the plasma ions, cx atoms of higher energies and with chaotically oriented velocities are generated. Although, namely, the latter are of our particular interest, to describe them firmly, all the species introduced have to be considered. In the vicinity of a circular opening in the wall (see Figure 1b), any neutral particle is characterized by the components V x and V r of its velocity, where x and r are distances from the wall and opening axis (see Figure 1b).

Recycling molecules and atoms
Recycling species, lunched from the wall, comprise backscattered atoms and desorbed molecules. Henceforth, we assume that the x components of their velocities have single values V bs, m and distinguish two groups of particles moving outward and toward the opening axis, with V r ¼ V bs, m and V r ¼ ÀV bs, m , correspondingly. The magnitudes of V bs, m are predefined by the generation mechanism for the species in question. Within the plasma the point with certain coordinates x and r can be reached only by particles, coming from the wall points 0; r þ x ð Þfor V r < 0 and 0; jr À xj ð Þfor V r > 0, respectively. By taking into account the attenuation processes with neutrals in the plasma, one gets n bs, m x; r ð Þ ¼ n bs, m 0; r þ x ð Þþn bs, m 0; where ν a ¼ n k a ion þ k a cx À Á and ν m ¼ n k m dis þ k m ion þ k m cx À Á are the decay frequencies for atoms and molecules, correspondingly, with n being the plasma density assumed the same for electrons and ions, k a, m ion are the rate coefficients of the ionization by electrons, k a, m cx those for the charge exchange with ions, and k m dis that for the dissociation of molecules.

Franck-Condon (fc) atoms
The fc atoms are born by the destruction of molecules and have positive and negative values both of V r and of V x . Henceforth, we distinguish four groups of these species with the velocity and m is the atom mass. The source density of fc atoms of each group is as follows: and depends on r through n m . The particle densities n fcAEAE , where the first subscript AE corresponds to the sign of V x and the second oneof V r , change along the characteristics AEx ∓ r ¼ const according to the continuity equation: where l is the length of the corresponding characteristics. The boundary conditions take into account that fc atoms are reflected with the probability R fc from the wall, x ¼ 0, and are absent far from the wall, x ! ∞. For the total density of fc atoms, n fc ¼ n fcþþ þ n fcþÀ þ n fcÀþ þ n fcÀÀ , one can obtain the Heaviside function Θ ξ ≥ 0 ð Þ¼1, Θ ξ < 0 ð Þ¼0 and r 0 being the opening radius (see Figure 1b).

Charge exchange cx atoms
The velocity distribution function of cx atoms, f cx x; r; V x ; V r À Á , is governed by the kinetic equation: Here, S cx ¼ S 0 cx þ S 1 cx is the total density of the source of cx atoms, with S 0 cx ¼ n k m cx n m þ k a cx Â n bs þ n fc À Á and S 1 cx ¼ nk a cx n cx being the contributions due to charge exchange collisions with ions of primary neutrals and cx atoms themselves, correspondingly, where is the density of the latter; the velocity distribution of the source is assumed as the Maxwellian one with the ion temperature T i , and is the ion thermal velocity.
2.3.1. Integral-differential equations for the density of cx atoms are governed by the equations, following from the integration of Eq. (3) over these ranges: with From the equation above, one gets the following ones for the variables The latter equation is straightforwardly integrated with respect to x, and for the r-component of the flux density, one gets Here, x 0 is the position, where the boundary condition is imposed; since there is no influx of cx atoms from the wall, γ l x ¼ 0 ð Þ¼0 for V x > 0, neutral species are absent deep enough in the plasma and γ l x ! ∞ ð Þ!0 for V x < 0. The found expression for γ l can be made more convenient if (i) the x-variation of terms involved is approximated by linear Taylor series: and (ii) by taking into account that the width of the region, occupied by cx atoms, is of several mean free path lengths (MFPL), that is, x À x 0 j j> V x j j=ν a typically. As a result we get By substituting the last expression into Eq. (5), this is reduced to the following one: Equation (7) can be integrated as the first-order equation with respect to x, with the boundary conditions similar to those for γ l , that is, one obtains the following integral-differential equation: where For the total density of cx atoms, one has n cx ¼ P lmax l¼1 η l .

Diffusion approximation
By neglecting the r derivative and by considering a single V r value, Eq. (8) is reduced to an integral one obtained in the 1D case (Eq. (16) in Ref. [6] or Eq. (11) in Ref. [4]). In this case l 1, δ l ¼ 1 and n cx ¼ η 1 . As it was demonstrated in Ref. [4], this integral equation can be reduced to an ordinary differential diffusion equation, if (i) the electron temperature is low enough and the ionization rate is much smaller than that of the charge exchange, k a ion ≪ k a cx , and (ii) the characteristic dimension for the plasma parameter change is much larger than the typical MFPL for cx atoms V th =ν a . The same procedure can be performed in the case under the consideration, providing the following 2D diffusion equation: with D th ¼ V 2 th = 2ν a ð Þ. The boundary conditions to Eq. (9) imply as follows: (i) cx atoms leave the plasma with the velocity close to the ion thermal one, and (ii) neutral species are absent far from the wall, that is, n cx x ! ∞ ð Þ¼0. Additionally, ∂ r η l ¼ 0 and ∂ r n cx ¼ 0 for r ¼ 0 and r ! ∞ for Eqs. (8) and (9), respectively.

Assessment of the velocity space integral I α
By inspecting Eq. (8), one can see the cause for large calculation time for kinetic computations for cx atoms. Even in the 1D 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 x ð Þ. In a 2D situation, the runs through the grid points in the r; V r À Á phase subspace are involved additionally into calculations. These can be, however, efficiently parallelized. The integral I α is appeared due to the generation of cx atoms with the ion Maxwellian velocity distribution and consequent attenuation by ionization and cx collisions. With a high accuracy, it can be assessed by an approximate pass method [12]. 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. The characteristic points u 1, 2 and u m are defined by the conditions F 0 α u m ð Þ ¼ 0 and F 00 α u 1, 2 ð Þ ¼ 0 for the derivatives F 0 α ¼ ÀF α g 1 =u 2 and F 00 For u m one can use the Cardano formula for the real root of a cubic equation By searching for u 1, 2 , we notice that for the roots of interest the equation g 2 u 1, 2 ð Þ ¼ 0 reduces to the following ones: With properly selected initial approximations accurate enough, solutions of these equations are found after three to five iterations with the Newton approach.
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 . The coefficients are selected to reproduce F αm;1; In addition, for u > u 2 the factor exp Àu 2 À Á leads to a very fast decay of F α with increasing u. Therefore, for the aim to assess the corresponding contribution to I α , we approximate F α by exp Àu 2 À α=u 2 À Á =u 2 . This results in the following expression: where δ 1, 2 ¼ u m À u 1, 2 j j . Figure 2b shows I α versus α found with the pass method, outlined above and evaluated numerically. By approximating the numerical results very accurately, the approximate pass method procedure allows to reduce the CPU time by a factor of 50 compared to direct estimates.
Alternatively to the usage of formula (10), one can calculate the integral I α in advance and save the results in a table. Then, for any particular α, the integral I α is interpolated from the values in this table. However, some time is needed to find the proper α interval, and this time grows up with α ! 0 since I α ! ∞ as Àlnα. Thus, the density of the data in the table has to be increased as α ! 0, and it is not obvious that such a way is more time-saving than formula (10).

Numerical procedure
The procedure to solve Eq. (9) numerically is organized as follows: For the problem in question, with the plasma profiles assumed unaffected by the processes in the vicinity of the duct opening, the integral I α x; y ð Þcan be assessed once in advance. Then, i. the source density S cx x; r ð Þis calculated, in the first approximation with n cx ¼ 0; ii. for all x the first term on the right-hand side of the equation is computed as a function of r, and the second-order ordinary Eq. (8) is solved by the method outlined in [10] for all V l ; iii. the next approximation for n cx x; r ð Þ and S cx x; r ð Þ is computed, and the procedure is repeated till the relative error in, for example, n cx 0; 0 ð Þ, becomes smaller than a desirable one.
A typical behavior of the error with the iteration number for calculations presented in the next section is shown in Figure 3. One can see that with consequent iterations the calculation error reduces steadily without any noise and the machine accuracy can be actually achieved.
Normally, a diffusion approximation is used to describe cx atoms to reduce CPU. The present approach to solve kinetic equation makes this unnecessary. To get the density of cx atoms, n cx x; r ð Þ, with a relative error below 10 À6 , kinetic calculations require much less CPU time than it is needed to solve the diffusion equation. Thus, for 500 and 30 grid points in the x and r directions, respectively, l max =10 and V max ¼ 3V th , the CPU times for a 1 GHz processor are 30 s and 112 min, correspondingly. Of course, the time for kinetic calculations increases significantly if neutrals affect the plasma background, and I α x; y ð Þ has to be evaluated by each iteration. Even in this case kinetic computations for cx atoms are done faster by factor of 4. Figure 4 shows the profiles of the plasma parameters, assumed homogeneous on the flux surfaces, across the surfaces, computed in [13] with the input parameters from the European DEMO project [1,2]. The distance x from the first wall to a certain flux surface depends on the poloidal angle θ (see Figure 1a) and in Figure 4 x corresponds to the torus outboard, θ ¼ 0, where this distance is minimal. The profiles of the cx atom density found with these parameters and for the duct radius r 0 ¼ 0:05 m, at the opening axis, r ¼ 0, and far from it, r ¼ 0:25 m are shown in Figure 5. There is a noticeable difference between the profiles of n cx calculated in the diffusion approximation (5a) and kinetically (5b). In particular, the difference between the density profiles at these two positions is significantly larger in the diffusion approximation than for those computed kinetically. In the latter case, n cx is still noticeable in much deeper and hotter plasma regions. This circumstance is of importance for the erosion of installations in ducts due to physical sputtering because this is very sensitive to the energy of impinging particles. To demonstrate this the erosion rate of a Mo mirror, imposed into the duct of three different radii r 0 , at the distance h from its opening (see Figure 1b) is assessed. To do such assessment, consider an infinitesimally thin toroid within   the plasma, with a square cross section of the width dr, thickness dx, and the radius r, situated at the distance x from the wall. Cx atoms with the energies within the range E, E þ dE are generated in the toroid with the rate:

Results of calculations
where is the Maxwellian energy distribution function of ions.
Since cx atoms move in all directions, some of them hit the mirror inside the duct (see Figure 1b). Henceforth, we analyze the surface position at the duct axis. One can see that cx atoms generated only in toroids with r ≤ r max ¼ r 0 1 þ x=h ð Þcan hit this mirror point. The contribution of the plasma volume in question to the density of the cx atom flux perpendicular to the mirror surface is is the cosine of the atom incidence angle with respect to the duct axis. The exponential factor in Eq. (13), with λ x; E ð Þ ¼ U x = ffiffiffiffiffiffiffiffiffiffiffi ffi 2E=m p , takes into account the destruction of cx atoms in ionization and charge exchange collisions on their way through the plasma. By ionization the atom disappears. By charge exchange a new cx atom is generated. The latter process is taken into account by the contribution S 1 cx in the source density S cx .
The energy spectrum of cx atoms, hitting the mirror, is characterized by their flux density γ ¼ Ð dj cx =dE in the energy range dE, where the integration is performed over r and x. By proceeding from r to s, according to the relation r ¼ The density of the outflow of mirror particles eroded by physical sputtering with cx atoms is as follows: Here, Y sp is the sputtering yield, whose dependence on E and s is calculated by applying semiempirical formulas from Ref. [14]. The erosion rate, measured henceforth in nm per full-power year nm=pfy ð Þ , is h sp ¼ Γ sp =n sp , with the particle density of molybdenum n sp . It is shown in Figure 6 versus the distance h between the wall position and the mirror point in question. One can see that in the former case the erosion rate is by a factor of 2 larger than in the latter one. It is of importance to notice that h sp , found with n cx computed kinetically, does not unavoidably excides that obtained with n cx calculated in the diffusion approximation. The results above have been gotten for a mirror positioned in a duct at the plasma outboard, that is, at the largest major radius R (see Figure 1a). Here, the local gradients of the plasma parameters are the largest because the distance between two particular flux surfaces has the minimum. Therefore, cx atoms penetrate, before they are ionized, into plasma regions with higher ion temperatures. The situation is different at the torus top where the corresponding distance has maximum and neutral species are attenuated already at significantly low n and T i . This leads to noticeably smaller h sp (see Figure 7). In this case, oppositely to the situation at the outboard, h sp is larger if n cx is computed in the diffusion approximation.
Due to technical requirements, an acceptable mirror erosion rate in a future fusion reactor should not exceed 1nm=pfy. As one can see in the case of a duct at the torus outboard, this target level is exceeded significantly for all h and r 0 under consideration. Seeding of the working gas into the duct is considered as a possible way to diminish the erosion rate below the maximum level allowed. The energy of cx atoms, coming into the duct, is reduced through elastic collisions with gas molecules, before cx atoms hit the mirror. Elastic collisions between neutral species lead to scattering on large angles. Consider a cx atom which, without gas in the duct, can hit the mirror directly. If the gas is seeded, in a narrow duct in question, with r 0 ≪ h, any collision of the cx atom with gas molecules leads to such a change of the atom velocity that with an overwhelming probability the atom will many times strike the duct wall before it gets the mirror. Through the collisions with the wall, the cx atom loses its energy so dramatically that it cannot contribute to mirror erosion. The reduction of the flux j E of incoming cx atoms with the energy E by collisions with the gas in the duct is governed by the equation: where σ el ≈ 3:8 Á 10 À19 E À0:14 m 2 ; eV Â Ã is the cross section for elastic collisions between cx atoms and molecules of hydrogen isotopes [15]. Consequently, the density of the outflow of mirror particles eroded by physical sputtering with cx atoms is modified, compared with Eq. (15), as follows: In Figure 8, the calculated dependences of the erosion rate h sp on h and r 0 are shown for several magnitudes of the density n g of the working deuterium gas in the mirror duct. One can see that the enhancement of n g above a level of 2 Á 10 19 m À3 should lead to the reduction of h sp below the target level of 1nm=pfy. The question, to what extent the local plasma parameters may be changed by the outflow of the gas from the duct, has to be investigated in the future on the basis of approaches developed in [16]. There, it has been demonstrated that the ionization of the gas outflowing into the SOL can lead to dramatic growth of the local density and cooling of the plasma to a temperature of 1 eV. Such cold dense plasma cloud can affect the transfer of cx atoms in the plasma near the opening in the wall. (Note that the h sp scale is different than in Figure 6).

Conclusions
The iteration approach to solve 1D kinetic equation for cx atoms, proposed decades ago [6], has been elaborated further to describe the transport of these species in a 2D geometry, in the vicinity of a circular opening in the wall of a fusion reactor. Unlike the Monte Carlo methods, this approach does not generate statistical noise so that calculation errors can be reduced to the level restricted by the machine accuracy. In order to perform calculations for a broad range of input parameters and do a thorough comparison with the results, obtained in the diffusion approximation for cx atoms [13], the solving procedure has been accelerated by a factor of 50, by applying an approximate pass method to assess integrals in the velocity space from functions, involving the Maxwellian velocity distribution of plasma ions.
The found possibility to speed up kinetic calculations is of importance, in particular, to perform firm assessments of the erosion rate of the first mirrors in future fusion reactors like DEMO. For a mirror located at the torus outboard, more accurate kinetic calculations predict by a factor of 2 higher erosion rate than the approximate diffusion approach. The erosion rate can be reduced very strongly either by putting the mirror duct at the torus top or by seeding the working gas into the duct. In the latter case, the elastic collisions with molecules in the gas reduce significantly the fraction of cx atoms which can hit and erode the mirror.
Nomenclature D l, th Diffusivity of cx atoms in kinetic and diffusion descriptions, correspondingly