Scrape-off layer modeling with kinetic or diffusion description of charge-exchange atoms

Hydrogen isotope atoms, generated by charge-exchange (c-x) of neutral particles recycling from the first wall of a fusion reactor, are described either kinetically or in a diffusion approximation. In a one-dimensional (1-D) geometry, kinetic calculations are accelerated enormously by applying an approximate pass method for the assessment of integrals in the velocity space. This permits to perform an exhaustive comparison of calculations done with both approaches. The diffusion approximation is deduced directly from the velocity distribution function of c-x atoms in the limit of charge-exchanges with ions occurring much more frequently than ionization by electrons. The profiles across the flux surfaces of the plasma parameters averaged along the main part of the scrape-off layer (SOL), beyond the X-point and divertor regions, are calculated from the one-dimensional equations where parallel flows of charged particles and energy towards the divertor are taken into account as additional loss terms. It is demo...


I. INTRODUCTION
By recombination of charged plasma components on the divertor target plates, limiters and vessel wall of fusion devices neutral particle species are generated.These recycle back into the plasma volume and participate in diverse interactions with electrons and ions, essentially affecting the plasma state. 1 By charge-exchange collisions with ions, the so-called c-x atoms of a broad velocity spectrum are produced so that a kinetic description of these particles seems to be mandatory.Such a description is performed usually by statistical Monte Carlo methods, 2 whose main problem is the high level of accident errors, "noise."This leads normally to a significant deterioration in the convergence of the plasma model coupled to the neutral module, and self-consistent neutralplasma simulations turn into an extremely time consuming numerical procedure.For example, the calculations of twodimensional plasma profiles in the scrape-off layer (SOL) and divertor, done in the framework of the ITER project, 3 required a central processing unit (CPU) time of months on the modern computers.Therefore, it is very attractive to use the reduced hydrodynamic approximations for neutrals which can be realized without "noise."Recently, 4 it has been proven, by describing the c-x atoms in a detached divertor, that the diffusion approximation, being formulated for the atom pressure, provides results very close to those obtained by integrating the kinetic equation.It is, however, not surprising for the conditions of very low temperatures of plasma components considered in Ref. 4. Indeed, in such a case, the rate coefficient for atom ionization by electrons, k ion , is much smaller than that of their charge-exchange with ions, k cx .Therefore, during the life time of an atom till ionization it undergoes a lot of charge-exchange collisions, leading to chaotic changes in the atom velocity, so that its motion is like diffusion.
Out of the divertor, in the main part of the SOL, the dominant source of neutral species is the plasma recycling on the first wall of the machine vessel.The temperatures of electrons and ions vary across the SOL very steeply, by approaching the level of hundreds of electron-volts near the separatrix.Under such conditions, a diffusion approximation for c-x atoms is not rigorously valid since k ion is comparable with k cx , and the characteristic scales for the plasma parameter variation in the radial direction-to the penetration depth of atoms.Therefore it is of a significant interest to assess the errors which are done, by applying approximate but fast diffusion models for neutrals in the SOL region.By keeping this aim in mind and in order to perform a broad parameter study, we apply in the present work a modeling approach "lighter" than the normally used.Presently, the 2-D kinetic Monte Carlo description for neutral particles and fluid transport equations for transfer of charged particles and heat in the plasma are normally applied to model axissymmetric tokamak devices, see, e.g., Ref. 5. Our approach relies on the averaging of original transport equations for the plasma components along the main SOL fraction.This procedure is substantiated by the results of 2-D calculations themselves, see, e.g., Ref. 6, showing the plasma parameter profiles rather flat along the flux surfaces in the SOL main fraction, out of the compact region near the X-point and divertor.Therefore they can be well represented by the values averaged along the main SOL.The variation of the averaged parameters across flux surfaces are governed by one-dimensional equations, deduced from the original 2-D ones, with the parallel particle and heat flows towards the divertor involved as additional loss terms.The kinetic modeling for c-x atoms is done by applying an iterative procedure elaborated in Ref. 7 to solve the equation for their velocity distribution function.Being formulated in the integral form, 7 this equation reveals clearly the origin of enormous CPU time spent on modeling of c-x atoms: double integrals over the normal and velocity spaces have to be calculated for every grid point at each iteration.It is demonstrated here that for the Maxwellian velocity distribution function of plasma ions, the integration over the velocity space can be done with a very high accuracy by an approximate pass method. 8This allows to reduce the CPU time at least by a factor of 30-50, compared to direct numerical computations of the integrals, and permits to perform an exhaustive comparison of the coupled neutral-plasma calculations done with either kinetic or diffusion descriptions of c-x atoms.However, also under such an extreme speeding up of computations, the modeling of c-x atoms remains the most time demanding part.Moreover, the relevancy of this procedure in 2-D and 3-D geometries is not demonstrated yet.This makes even more important the validation of the diffusion approximation, which allows reducing the CPU time by orders of magnitude further and can be easily generalized on multi-dimensional situations.
The remainder of the paper is outlined as follows.In Section II the set of equations used to model the neutral species generated by the plasma recycling at the vessel wall is formulated.The pass method to calculate integrals over the velocity space, being involved in the kinetic equation in the integral form and the reduction in the limit k ion ( k cx of the kinetic description to the diffusion one, formulated for the c-x atom pressure, are also demonstrated here.The averaging of plasma transport equations along the flux surfaces in the SOL main fraction and the assessment of the loss terms due to particle and heat flows towards the divertor are performed in Section III.In Section IV computations are done for the input parameters from the European DEMO project 9,10 and the results obtained with either kinetic or diffusion description of c-x neutrals are compared.It is demonstrated that practically for all characteristics of interest, the diffusion approximation is sufficiently satisfactory.Conclusions are drawn in Section V.

A. Species included for consideration
The span of neutral species of hydrogen isotopes released from the first wall into the SOL plasma is defined by the phenomena at the wall surface.Electrons and ions leave the plasma due to transport processes and recombine into atoms on the wall.With the probability R i bs , being a function of the ion energy E i , and the mass and charge of the wall particles, these atoms escape back into the plasma, with an energy E i bs ՇE i .The rest fraction, 1 À R i bs , of atoms give up their energy to the wall particles and recombine into molecules.Finally, the molecules of hydrogen isotopes can be released from the wall with the probability x rec .A similar situation takes place for atoms generated in the plasma by chargeexchange collisions and escaping onto the wall; R cx bs and E cx bs denote the back-scattering probability and the energy of these species.The atom species above vanish in the plasma because of the ionization and charge-exchange collisions with electrons and ions.These are characterized, respectively, with the rate coefficients k a ion , being a function of the electron temperature, and k a cx , weakly depending on the ion temperature and atom energy. 11The decay of the densities of the backscattered atom species, n i bs and n cx bs , correspondingly, with the distance x from the wall is described by the relations where C p is the density of the charged particle outflow from the plasma onto the wall, j cx ðxÞ the flux density of c-x atoms; V i;cx bs % ffiffiffiffiffiffiffiffiffiffiffiffiffiffi E i;cx bs =m q are the velocities of the back-scattered atom components with m being the atom mass; i;cx bs ¼ ðk a ion þk bsÀi;cx cx Þn are the total frequencies of the atom species destruction, with the subscripts ðbs À iÞ and ðbs À cxÞ meaning k a cx calculated for the energy of the corresponding atom species, see the Appendix; n is the plasma density assumed the same for electrons and ions.
Neutral molecules released into the plasma are disintegrated in processes of dissociation and ionization by electrons and charge-exchange with ions, characterized by the rate coefficients k m dis ; k m ion and k m cx , respectively.For the molecule density n m one has and n À f c , respectively.The change of n 6 f c with x is governed by continuity equations is the characteristic velocity of f-c atoms.The difference of these equations for n þ f c and n À f c , respectively, provides where is the total density of f-c atoms.By differentiating the latter relation once more and substituting into the sum of equations above for n þ f c and n À f c , one gets a "diffusion" equation for n fc The boundary condition for this equation at x ¼ 0 follows from the back-scattering condition at the wall, and @ x n f c ¼ 0 at the plasma axis, x ¼ r w , where r w is the wall minor radius, i.e., the distance from the plasma axis to the wall.The distribution function of atoms generated by chargeexchange collisions with ions, u cx , is a function of x, the xcomponent of the atom velocity, V x , and the kinetic energy E ? of the atom motion in the plane perpendicular to x.The integration of the kinetic equation for c-x atoms with respect to E ?results in the following 1-D equation for the distribution function where cx ¼ nðk a ion þ k cx cx Þ, S cx is the density of the total source for c-x atoms with being the source part due to charge-exchange of primary neutrals, i.e., molecules, f-c atoms and species backscattered from the wall is the density of c-x atoms.The adequacy of the Maxwellian velocity distribution for the charge-exchange source term in Eq. ( 5), with the ion thermal velocity V th ¼ ffiffiffiffiffiffiffiffiffiffiffiffi ffi 2T i =m p , has been discussed in Ref. 7. Arguments, supporting this approximation, are given also in the Appendix of the present paper.
By considering separately the c-x atoms with V x 90, one can find the integral relations for the corresponding parts of the velocity distribution function and where By integrating f cx with respect to V x , one gets an integral equation for n cx , where the integral with and a xy ¼ jU x À U y j=V th ðyÞ.By taking into account the balance of hydrogen particles in the wall, one can find the density of the molecule influx where the outflows of f-c and c-x atoms are given by the relations

B. Assessment of kinetic integrals
An inspection of Equations ( 10)-( 12) allows to see the origin of large CPU time expenditures on kinetic calculations for c-x atoms: for each grid point x, one has to calculate the 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Þ and the plasma parameters.However, by looking at the shape of the function under the integral displayed in Figure 1(a) together with functions f 0;À1 ðuÞ for a xy ¼ 1, one can notice: in the case under consideration, the usage of an approximate pass method 8 can be sufficiently precise and bring large time savings compared to direct numerical calculations of I k .For each of the functions f k ðuÞ, one can distinguish 4u-intervals where it behaves principally differently: and u 2 < u, see Figure 1(a) where the points u 1;2 and u m are marked for k ¼ 1.These points are defined by the conditions For k ¼ 0, we get u m ¼ ða=2Þ 1=3 ; for k ¼ 1 and k ¼ À1 u m is a root of a cubic equation.In the former case, this equation has a single real root . In the latter one the expression for u m above is applicable but with b ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi a 2 À 1=216 p and for a !ffiffiffiffiffiffiffiffiffiffi 2=27 p only.A more complex analytical expression exists also for a < ffiffiffiffiffiffiffiffiffiffi 2=27 p .However, it is simpler and faster to use the Newton method 12 to find an approximate solution for g 1 ¼ 0. By taking as the initial approximation 1 ðu m Þ to find u m with a relative accuracy 10 À4 for any a.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: 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 k , we approximate the functions f k ð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 and where dx being the exponential integral function.Figure 1(b) shows I 1;0;À1 versus a 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 (14), one can calculate in advance the integrals I k ðaÞ as functions of a and save this in a table.Then, by solving Equation (11) for any particular a, the integral I k ðaÞ can be interpolated from the values in the table for arguments larger and smaller than that in question.However, the searching of the correct a-interval in this table requires also some time.In the case of the most important I 1, this time grows up with a !0 since I 1 ðaÞ ! 1 as E 1 ðaÞ 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 the formula (14).

C. Diffusion approximation
With the usage of the pass method for the assessment of the integrals I k , one-dimensional kinetic calculations for neutrals stand alone require 6-7 iterations to get a solution with an accuracy higher than 1%.It requires less than 3 s of the CPU time on a 1 GHz processor compared to 2 min needed for the same problem, being solved either by calculating the kinetic integrals numerically or with a Monte Carlo approach. 13For the accuracy 10 À6 adopted in the present study, kinetic neutral calculations require roughly 60 iterations.A converged solution of the coupled neutral-plasma FIG.
1.The functions f k ðuÞ ¼ exp ðÀu 2 À a=uÞu Àk for a ¼ 1 and k ¼ 1; 0; À1; the points u m and u 1;2 defined by the conditions f 0 k ðu m Þ ¼ 0 and f 00 k ðu 1;2 Þ ¼ 0, correspondingly (here for k ¼ 1) (a); the "kinetic" integrals I k ðaÞ ¼ Ð 1 0 f k ðuÞdu estimated by the pass method (solid curves) and computed numerically (dashed curves) (b).problem needs a reasonable CPU time of 30 min.Nonetheless, it is unclear whether it is possible to utilize a similar procedure in a 2-or 3-D geometry.Normally in such a case, a diffusion approximation is used for c-x atoms to avoid a huge noise.To derive this, the velocity distribution function of the particles in question is assumed as a shifted Maxwellian one and the fluid like equations are deduced for the continuity of particle density and momentum.Henceforth, by neglecting the mass velocity contribution to the inertia term, a diffusion like relation for the flux can be obtained, see, e.g., Refs.4, 14, and 15.However, the validity of the assumptions and steps described above are not strictly proved.Therefore, we deduce here the particle flux density in a diffusion approximation directly from the velocity distribution function for c-x atoms given by the relations ( 8) and (9).By taking into account that x and U x are uniquely related to each other, we introduce the variable As per definition, a diffusion approximation is applicable if the mean free path length (mfpl) of particles in question is much smaller than the characteristic dimensions for the change of the parameters involved.The typical mfpl for c-x atoms is k cx ¼ V th =ðk cx cx nÞ, but their total penetration depth is defined by k a ion .Thus, the diffusion approximation is relevant if k a ion ( k cx cx .Moreover, it is required (i) jV x j ( U x;y for typical V x ; x; y and (ii) and the Taylor's expansion Þis applicable.By substituting the expansion above into the relations (8) and ( 9), one can easily perform the integration over U y .By neglecting exponentially small terms proportional to exp ðÀU x =jV x jÞ, we get Because G is an even function of V x , the former term in f cx ðx; V x Þ does not contribute to the flux density of c-x atoms and we obtain By integrating the kinetic equation ( 5) with respect to V x , we get S cx0 þ nk cx cx n cx À nðk a ion þ k cx cx Þn cx on the right hand side.Thus the term nk cx cx n cx is cancelled out, and one has a continuity equation for c-x species as follows: ion nn cx : In the diffusion limit, the c-x atom mfpl k cx has to be much smaller than the characteristics change scale, in particular, for S cx0 .Therefore the contribution from the second term in j cx , being of k 2 cx @ 2 x S cx0 , can be neglected compared to the term S cx0 on the right hand side.Thus, we get the following diffusion equation for the density of c-x atoms The boundary condition to Equation ( 16) at the wall implies that c-x atoms leave the plasma with the thermal velocity, i.e., j cx ð0Þ ¼ Àn cx r ffiffiffiffiffiffiffiffiffiffi T i =m p at x ¼ 0, with the factor r % 1:07 seems to be the best choice; j cx ¼ 0 at the plasma axis, x ¼ r w .

A. Basic transport equations
In the SOL, the charged plasma components, electrons, and fuel ions, are described by two-dimensional fluid equations.These govern the variations of the plasma density n, electron and ion temperatures T e and T i , correspondingly, across the flux surfaces, the coordinate x, and along the magnetic field, the direction l, where C ?;jj are components of the particle flux density along the x and l directions, q e;i ?¼ Àj e;i ?@ x T e;i þ 1:5C ?T e;i and q e;i jj ¼ Àj e;i jj @ l T e;i þ 2:5C jj T e;i the components of the heat flux densities, with j e;i ?;jj being the components of the heat conduction; , the contributions to the plasma particle source density due to ionization of molecules and atoms, respectively; E m ion and E a ion the energy spent by electrons in ionization collisions; is energy transported to the ions by ionization and charge-exchange of atoms, with being the heat capacity density of c-x atom; the last terms on the right hand sides of the heat balance equations, Q ei , are due to the heat exchange between electrons and ions through coulomb collisions.
The boundary conditions to Equations ( 17)-( 19) at the wall, x ¼ 0ðr ¼ r w Þ, are fixed e-folding lengths of the parameters in question, d n ¼ @x=@lnn, d e;i ¼ @x=@lnT e;i ; at the separatrix, x ¼ x s ðr ¼ r s Þ, we prescribe the plasma density nðr s Þ and the densities of the heat outflows from the confined plasma region, q e;i ?ðr s Þ.Some fraction of neutrals, recycling from the wall, may penetrate through the SOL region into the confined plasma, 0 r r s .To take these into account, the density and temperature profiles have to be postulated there.We distinguish the edge transport barrier, r b ¼ r s À D b r r s , where with k n ¼ ln½nð0Þ=nðr b Þ and k e;i ¼ ln½T e;i ð0Þ=T e;i ðr b Þ.The parameters r b , r s , nð0Þ; nðr b Þ, nðr s Þ; T e;i ð0Þ; T e;i ðr b Þ are prescribed according to the European DEMO project, 9,10 T e;i ðr s Þ are got by solving the transport equations in the SOL.

B. Equations for parameters averaged along the SOL
As it was already stated above, we consider here the main part of the SOL, henceforth the main SOL, out of the divertor and the X-point vicinity where a significantly enhanced interaction between the plasma and neutrals often takes place.A good estimate for the extent of the main SOL along the magnetic field is two connection lengths L between the high and low field sides of the torus.To obtain equations, describing the variation in the x-direction of the parameters averaged along the main SOL, we integrate Equations ( 17)-( 19) with respect to l in the range 0 l L, with l ¼ 0 corresponding to the symmetry plane where C jj ; q e;i jj ¼ 0. As a result one gets dhq e ?i=dx ¼ ÀE m ion hS m ion i À E a ion hS a ion i À hQ ei i À q e jjL =L; (21) where hÁ Á Ái ¼ Ð L 0 Á Á Á dl=L, C jjL and q e;i jjL are the densities of the parallel losses of charged particles and their energy at the end of the main SOL.
To relate C jjL and q e;i jjL to the SOL averaged parameters, hni and hT e;i i, we apply an approximate method from Ref. 16  for solving the parabolic partial differential equations, e.g., a diffusion one.This method is based on the concept that equations in question allow neither the periodic sign-changing solutions nor a similar behavior of individual terms in them.It presumes that the solution profile in a certain direction is mostly controlled by the corresponding transport term and is not very sensitive to the spatial variation of other terms.By following this presumption, we rewrite Equations ( 17)- (19)  as follows: @ l C jj ¼ S; @ l q e;i jj ¼ W e;i ; with some unknown source/sink terms S and W e;i .By assuming these independent of l and integrating equations above from l ¼ 0, we get C jj ðlÞ=C jjL ¼ q e;i jj ðlÞ=q e;i jjL ¼ l=L: With the relations for q e;i jj above and by adopting the Spitzer-H€ arm approximation for the parallel heat conduction, 17 j e;i jj % A e;i jj T 2:5 e;i , with A e jj % 10 18 m À1 s À1 eV À2:5 and A i jj % 3:3 Â 10 16 m À1 s À1 eV À2:5 , one obtains for the temperatures T 2:5 e;i dT e;i 1 À T e;i =T c e;i ¼ q e;i jjL dl 2 2LA e;i jj : Here T c e;i ¼ q e;i jjL =ð2:5C jjL Þ are the maximum temperature values, corresponding to the whole heat flux transported with convection by the plasma particle flow only.The latter equation can be analytically integrated, providing a transcendental algebraic equation for T e;i ðx; lÞ, see Ref. 18. Here, however, we prefer to represent the result of integration as a polynomial algebraic equation for h ¼ T e;i ðl; xÞ=T e;i ð0; xÞ, in the following form convenient for a fast finding of the root in question with the Newton method: The parameters h L ¼ T e;i ðL; xÞ=T e;i ð0; xÞ and n 0 ¼ T e;i ð0; xÞ= T c e;i ðxÞ are external for the problem under consideration; these can be defined by including the X-point and divertor regions into the analysis.However, for the conditions of either strong recycling or detachment, expected in the DEMO divertor, heat will be transported in the SOL mostly by the heat conduction.Thus n 0 , h L ( 1 and the summations in the equation above converge fast.Figure 2(a) shows the average value of h along the main SOL, hhi, versus the parameters h L and n 0 ; it is in a very narrow range of 0:86 À 0:87.For the relation between the heat loss density and the temperatures of the plasma components average along the main SOL, hT e;i iðxÞ ¼ T e;i ð0; xÞhhi, we get q e;i jjL ¼ vA e;i jj hT e;i i 3:5 =L; (24 The dependence v on the parameters h L and n 0 is displayed in Figure 2(b); roughly v % 1 þ 2:5n 2 0 for small n 0 .In concluding this section, we assess the relation of the plasma particle loss density C jjL to the plasma parameters averaged over the main SOL.To maintain a steady state, the outflow from the SOL into the divertor has to be exhausted, by pumping out a part x pump ( 1 of neutrals born by plasma recombination at the target plates, i.e., C jjL ¼ x pump C jjp , with the plasma flux to the plate 6 where n p and M p ! 1 are the plasma density and Mach number, correspondingly, at the target; due to high frequency of coulomb collisions, electrons and ions have here roughly the same temperature T p .In the main SOL, the parallel plasma momentum is lost with the perpendicular transport and friction by charge-exchange collisions with recycling neutrals.All together, these losses are characterized by the time s loss $ 10 À4 s and the motion equation is as follows: jj =n being the total plasma pressure.By integrating Equation (25) along the main SOL, 0 l L, with relations (23) taken into account, one gets where Pð0Þ ¼ nð0Þ½T e ð0Þ þ T i ð0Þ and PðLÞ are the total pressure values at the SOL symmetry plane and at the entrance of the divertor region, respectively.P is normally reduced further in the divertor by a factor d d < 1, the so-called "detachment degree," and at the target it approaches the level By introducing g ¼ n=nð0Þ and combining Equations ( 26) and ( 27), we get and the "divertor impact factor" To make a guess about the latter, we apply the parameters foreseen for ITER, 3 T p % 3 eV, M p ¼ 1, d d ¼ 1=2, and x pump % 10 À2 .This results in V d % 1:5 Â 10 6 ms À1 , exceeding significantly the ion sound velocity of 10 5 ms À1 for T e;i % 200 eV typical of the DEMO SOL, see Figure 3. Since l; hgihhi $ 1, the flow in the main SOL part is deeply subsonic, as it was assumed and employed above, by the assessment of the parallel heat transport.Thus, v ¼ 1 can be used in the relation (24) for q e;i jjL .To assess hgi explicitly, we integrate Equation (25) from 0 to l and obtain or a quadratic equation for gðlÞ, with the solution of interest: Taylor's expansion of the square root term can be done and, finally, we get hgi % We have to admit that the divertor impact factor V d is probably the most obscure among those involved in Equations ( 28) and (29).Therefore the influence of V d will be examined by varying it in a broad range.

IV. RESULTS OF CALCULATIONS
Below we present the results of calculations for the "standard" set of the input parameters taken from the European DEMO project: 9,10  ; T e ð0Þ ¼ 27:4 keV; T e ðr b Þ ¼ 8 keV and T i ¼ 1:5T e is assumed in the plasma core, 0 r r b ; the heat flux densities at the separatrix q e;i ?ðr s Þ ¼ 54 kW m À2 ; the connection length in the SOL L jj ¼ 90 m, the divertor impact factor V d ¼ 1:5 Â 10 6 ms À1 and the e-folding lengths of the plasma parameters at the wall of tungsten, d n;e;i ¼ 10 À3 m.In the SOL, we assume a diffusive perpendicular plasma particle transport, C ? ¼ ÀD ?@ x n, with the Bohm's diffusion and heat conduction D ?¼ D B T e =ð16eB T Þ; j e;i ?¼ 1:5nD ?; at the magnetic field of 5:7 T. The rate coefficients for inelastic collisions of neutral species with electrons and ions, k m dis ; k m;a ion and k m;a cx , are taken from Ref. 11. Figure 3 shows the radial profiles of the SOL averaged values for the particle densities of molecules, hn m i, all atom species, hn a i, and plasma, hni; the densities of sources for cx atoms, hS cx i, and charged particles, hS ion i; the temperatures of electrons, hT e i, and ions, hT i i.The results of calculations with the kinetic and diffusion descriptions of c-x atoms are displayed by dashed and solid curves, respectively.The profiles of hn m;a i and hS cx;ion i, calculated with two approaches, differ maximally by 15%.The difference in the plasma parameter profiles is much smaller and can be neglected in view of uncertainties, e.g., in the perpendicular transport coefficients.Indeed, by varying the parameter a B ¼ D ?=D B from 0.8 to 1.2, we get a significantly larger change in the profiles, as it is demonstrated in Figure 4 for hT e;i iðxÞ.As it was mentioned above, the divertor impact factor V d , characterizing the losses of charged particles along magnetic field, is probably the most vulnerable input quantity for the one-dimensional Equations ( 20)-(24).In Figure 5, we demonstrate the profiles of the plasma component temperatures for different V d values calculated with either kinetic or diffusion description of c-x atoms.One can see that the decrease of V d leads to a noticeable drop of the temperatures close to the separatrix due to the enhancement of the heat loss with  the parallel convection of plasma particles.Although these changes are not very impressive, a firmer modeling approach has to include a coupling with a 2-D model for the X-point and divertor region, predicting reliably V d as a function of the x-coordinate.The relatively moderate effect of the V d value on the hT e;i iðxÞ profiles allows to anticipate that an iterative coupling procedure can work successfully.
For the life time of the reactor wall, the energetic spectrum of c-x atoms, leaving the plasma, is of importance because such particles can have energies significantly exceeding the temperatures of the plasma components near the wall.This spectrum we characterize by the density c cx ðEÞ of the outflow of c-x atoms with the total energy E. It can be calculated from the distribution function u cx ðx; . Exactly the same expression follows by calculating c cx ðEÞ for a given source density of c-x atoms and the Maxwellian energy distribution function for ions, an approach which can also be applied if the density of c-x atoms is calculated in the diffusion approximation. 19The energy spectrum of the atom outflow computed with both kinetic and diffusion description of c-x atoms for three magnitudes of the divertor impact factor V d is demonstrated in Figure 6.One can see that the diffusion approximation systematically overestimates c cx ðEÞ in the vicinity of its maximum but underestimates for higher energies.The latter is also mimicked in the heat load density on the wall due to the outflow of c-x atoms In Figure 7 q cx is displayed as a function of V d calculated for different magnitudes of the parameter a B ¼ D ?=D B , with FIG. 6.Energy spectra of the outflow of c-x atoms to the wall for different magnitudes of the divertor impact factor, calculated with the kinetic (dashed curves) and diffusion (solid curves) descriptions of c-x atoms.either kinetic or diffusion description of c-x neutrals.The underestimation in the diffusion approximation of c cx ðEÞ for large E results in q cx smaller by maximally 20% À 30%.This error can be significantly reduced if after the convergence of coupled plasma-neutral calculations, the final computation of c-x atoms is done kinetically.The dependencies of q cx on V d and a B found in this case are also shown in Figure 7.

V. CONCLUSIONS
The main aim of the present study is to validate the diffusion approximation to describe atoms generated by the charge-exchange of neutral species recycling from the wall of a fusion reactor.The demand to apply such an approximation is dictated by the high level of accident errors, "noise," in the characteristics of neutral species computed by the statistical Monte Carlo methods, leading to a slow convergence of coupled neutral-plasma computations.The present consideration is confined to the main SOL, and equations for the variation of the parameters with the distance from the wall are integrated.It is demonstrated that the solving of the kinetic equation for c-x atoms, represented in the integral form, can be accelerated at least by a factor of 30 if the arising integrals over the velocity space are evaluated by an approximate pass method.A possibility to apply this approach to 2-D or 3-D situations has not been proven yet.Therefore, the credibility of the diffusion approximation, being strictly valid for low temperatures when the ionization rate coefficient k a ion is much smaller than that of the chargeexchange, k a cx , has been thoroughly investigated for the SOL conditions where k a ion $ k a cx .By applying this approximation, which can be easily generalized for 2D and 3D configurations, calculations for c-x atoms can be speeded up additionally by a factor of 10 3 .To perform a self-consistent modeling of the plasma in the SOL, equations for the variation across the flux surfaces of the density and temperatures of electrons and ions, averaged along the SOL, are deduced by integrating the two-dimensional fluid equations along the magnetic field.It is demonstrated that the parallel outflows of heat from the integration domain towards the divertor can be firmly assessed and related to the parameters averaged along the main SOL only.In the case of plasma particle losses, the so-called "divertor impact factor" has to be introduced into the model additionally.
Coupled 1D models for neutral and charged species, either with kinetic or with diffusion description of c-x atoms, have been numerically realized and calculations performed for the input parameters from the European DEMO project.It is shown that the diffusion approximation works extremely good, by providing practically the same profiles across flux surfaces for the SOL averaged plasma parameters as those obtained with the kinetic description for c-x atoms.The maximum discrepancy between the two approaches appears for the characteristics of c-x species themselves.So their energy flux onto the wall is underestimated with the diffusion approximation by 20%-30%.Nonetheless, this discrepancy can be significantly reduced if the plasma parameters found in the latter case are used for the final kinetic calculation of c-x atoms.Our investigation permits to validate firmly the diffusion approximation for c-x atoms for the SOL conditions of relatively high temperatures with comparable rates of ionization and charge-exchange and allows an enormous saving of the CPU time by several orders of magnitude.To our mind, this conclusion can also be practical for 2-D and 3-D geometries; however, careful peculiar considerations have to be performed also in these situations.
In the present study, an idealized situation with a clearance between the wall and separatrix constant along the SOL has been assumed.In a realistic case, this changes along the flux surfaces due to shift, elongation, and triangularity of the surfaces.Nonetheless, we want to stress that the onedimensional models for neutral and plasma parameters in the main SOL are also relevant to these conditions.The penetration depths for neutral species are normally much shorter than the characteristic distance where the magnetic geometry changes in the poloidal direction.Thus, a one-dimensional description of neutrals, recycling from the wall, can be applied individually at different poloidal positions, see Ref. 20.In Equations ( 20)-( 22) for charged particles, the terms S m;a ion and S cx , arising due to the presence of neutrals, are linearly related to their densities.If the latter are calculated in the one-dimensional approximation for poloidal positions with noticeably different distances between the flux surfaces, one can straightforwardly perform averaging in the equations for the plasma parameters, if these change weakly along the main SOL.

APPENDIX: VELOCITY DISTRIBUTION OF c-x ATOM SOURCE
In the kinetic Boltzmann equation for c-x atoms, 7 the source term due to collisions with ions is as follows: where v and v 0 are the velocities, u i and u a , the velocity distribution functions of ions and atoms, respectively.In a very broad range of the energy E ¼ mjv À v 0 j 2 =4 of ion-atom relative motion, 1 eVՇEՇ10 keV, the rate coefficient of atom charge-exchange collisions changes very weakly, k a cx ðEÞ $ E 0:3 .Thus, without a large error, one can use for this the values, computed for the characteristic energies of interacting ions and atoms, E i ¼ 1:5T i and E a , correspondingly, with the latter being equal to E i bs and E cx bs for backscattered species, E fc and E cx ¼ W cx =n cx for f-c and c-x atoms where W cx is the thermal capacity of c-x atoms defined after Equation (19).With an accuracy of 2%, where k a cx is measured in m 3 s À1 , E i , E a -in eV and A i is the atomic weight of the hydrogen isotope in question.
With this approximation, we get The cross-section of charge-exchange is much larger than that of the elastic collisions between hydrogen ions and atoms.Therefore, the ion velocity distribution function is practically unperturbed by the charge transfer.By choosing u i ðvÞ, we take into account that (i) the ion thermalization time in the DEMO SOL will be of 0:1 ms, i.e., much shorter than the characteristic transport times along and across the SOL of 10 ms, see Section III B; (ii) ions are magnetized so their mass velocity in the x-direction across the magnetic field is much smaller than V th .Thus the Maxwellian distribution function can be assumed for ions By integrating C cx over E ?, one gets the source term assumed in Equation ( 5), S cx .

5 ;
(3)   where j m ð0Þ is the molecule influx density and V m is their mean velocity.In collisions of molecules with the plasma particles, Franck-Condon (f-c) atoms, with a characteristic energy E f c ' 3:5 eV, are generated.The velocity of f-c atoms is chaotically directed, and one can introduce separately the densities of particles moving from and towards the wall, n þ f c Phys.Plasmas 23, 122512 (2016) e b nðr b Þ; T e;i ¼ e s T e;i ðr s Þ þ e b T e;i ðr b Þ; with e s ¼ ðr À r b Þ=D b ; e b ¼ 1 À e s , and the core, 0 r r b , where n ¼ nð0Þ exp ðÀk n r 2 =r 2 b Þ; T e;i ¼ T e;i ð0Þ exp ðÀk e;i r 2 =r 2 b Þ; FIG.2.The SOL averaged value hhi of dimensionless temperature (a) and the factor v in the relation (23) (b) versus h L calculated for different magnitudes of the convection fraction: n 0 ¼ 0:05 (solid line), 0.1 (dashed line), 0.2 (dotted line), and 0.3 (dashed-dotted line).

FIG. 3 .
FIG.3.The radial profiles of the SOL averaged molecule and atom densities (a), the source densities of c-x atoms and charged particles (b), plasma density (c) and electron and ion temperatures (d), calculated with the kinetic (dashed curves) and diffusion (solid curves) descriptions of c-x atoms.

FIG. 5 .
FIG.5.The radial profiles of the SOL averaged temperatures of electrons and ions calculated for 3V d and 0:3V d with kinetic (dashed curves) and diffusion (solid curves) descriptions of c-x atoms.