Study of the erosion and redeposition of W considering the kinetic energy distribution of incident ions through a semi-analytical model

In today’s nuclear fusion devices, erosion of high-Z metallic plasma-facing materials (PFMs) is mainly caused by physical sputtering. That is, by the exchange of energy between plasma ions and the atoms in the walls. In most of the numerical codes currently in use impinging plasma is approximated as a fluid. By averaging the incident particles’ energy distribution the high-energy population of the eroded material is underestimated. For heavy materials such as W, high-energy eroded particles tend to ionize far from the wall and they are less affected by the sheath electric field hence, not being attracted back to the wall, they have a higher chance to contaminate the core plasma. This could in turn result in an underestimation of the net erosion sources. In this work, a semi-analytical model was developed to include the energy distribution of the incident particles. Then, by Monte Carlo method, the net erosion of tungsten from a smooth PFM was calculated. The results show that the kinetic description in energy is important only for incident particles ionized once. For instance, it is particularly important for plasma ions such as Deuterium. It is seen that Deuterium contribution to the W net sources is not always negligible if compared to light impurities or to tungsten self-sputtering in the range of plasma parameters tested. Finally, results show that the difference between the fluid and kinetic models becomes more pronounced for high-screening plasma conditions.


Introduction
In the future, magnetic nuclear fusion reactors might rely on metallic plasma-facing materials (PFMs) such as tungsten (W), due to their low sputtering yield, high melting point, and low tritium retention [1]. In these reactors, accurate assessment of the impurity source is fundamental to the operation of high-Z plasma PFMs, as the penetration of heavy impurities into the core leads to unacceptable radiative losses due to their high cooling factor [2]. This also affects plasma transport and stability thus, the overall performance of the tokamak [3]. Different plasma-wall interaction (PWI) mechanisms underlie wall erosion. For PFMs composed of heavy atoms, the predominant mechanism is the so-called physical sputtering [4]. The total flux of eroded particles is called the 'gross erosion flux' (Γ gross ), which can be expressed as the incident flow (Γ inc ) times an effective sputtering yield (Ȳ). This quantity is experimentally measured via spectroscopy measurements [5]. As shown by [6], not all gross flux is a source of impurities for the plasma. In fact, a fraction of the sputtered particles readily redeposits on PFMs as a consequence of their gyromotion. In [6] the fraction of eroded and redeposited particles was defined as the prompt fraction f p (non-prompt fraction f np = 1 − f p ), and it was used to adjust the sputtering yield redefined as 'net' sputtering yield. This led to a distinction between the gross flux Γ gross = Γ inc ×Ȳ and the net flux Γ net = Γ gross × f np . It was shown by [7] that the definition of f np adopted in [6] was not appropriate for high-Z PFMs. In particular, through a set of local Monte Carlo simulations, they noted that f np was largely overestimated if the effect of the plasma sheath electric field was not considered as done in [6]. For this reason in [7] f np was redefined as f non−redep , that is, the fraction of all non-redeposited particles within the Monte Carlo simulation. The importance of the electric field in the plasma sheath was reaffirmed in [8], where its predominance was confirmed even in the presence of thermal forces. Moreover, in [9], it was noticed that the net flux of particles is mostly populated by the fast tail of eroded particles. Hence, the resulting f non−redep , was affected by the angular and energy distribution of sputtered particles. This is because of essentially two reasons: (i) faster neutral particles tend to ionize further from the wall, thus, in a zone of low electric field strength, (ii) faster particles have a higher chance to overcome the electric field potential in case they are ionized inside the sheath. So it is clear that in order to describe the net sources a good description of the sputtered particles' energy and angular distribution is needed. In current codes such as ERO [10], ERO2 [11], DIVIMP [12], IMPMC [13], and GITR [14] the sputtered particle flux distribution is generally represented with a Thompson distribution [15] or by some variants adopted for the impact of low-energy plasma particles ( [16,17]). These distributions are calculated experimentally through mono-energetic beams of incident particles. While they are applied, impinging ions are approximated with a mono energetic beam whose energy E i is expressed as a function of integrated or 'fluid' plasma quantities, e.g. E i = 2T i + ∆ϕ Z, with ∆ϕ being the potential drop in the plasma sheath, see [18]. Neglecting the kinetic distribution by applying averaged 'fluid' quantities has its limitations, for example, in [19] and [14] it was shown that considering kinetic effects leads to a refined estimate of the sputtering yield for gross erosion. In this work, the goal is to go a step further by including the kinetic energy distribution of impinging ions in the assessment of net erosion. Results show that accounting for the whole energy distribution of impinging ions on W is mostly relevant when the incidence energy is close to the sputtering threshold energy limit. In general, only impurities ionized once and plasma ions were close to the sputtering threshold in the simulation here tested. Comparisons between the W net source related to once-ionized impurities and Deuterium (D), the contribution of D is not always negligible. Moreover, it is discussed the dependence of f non−redep with respect to the magnetic field vector and the applied electric potential drop model inside the sheath. Finally, the energy distribution of incident particles shows a good agreement with the one obtained via Monte Carlo simulations.

Sheath model
A Deuterium plasma background is considered, plasma oscillation modes and turbulent fluxes are neglected if compared to convective particle flux. It is considered a steady state flux of impinging ions distributed as f i (⃗ v) at the sheath entrance (SE) accelerated inside the plasma sheath. Ions are immersed in a uniform magnetic field B oblique with angle α B respect to an ideal flat surface (no ∇ ⃗ B forces and roughness considered) see figure 1. The electric field is considered one-dimensional and normal to the surface. Ions at the SE respect the Bohm criterion and their energy gain both inside the scrape-off layer (SOL) and the sheath is supposed to follow the rule proposed in [20]. Reflections of impinging ions on the target material and redeposition-related sputtering are neglected. The sheath is supposed collisionless and isothermal T e ≈ constant, with no gradients along x, y directions. These are working assumptions, in reality, because of charge-exchange collisions the sheath presents a certain degree of collisionality especially for low electron temperature at the sheath edge (T e < 10 eV), and densities n e > 10 19 m −3 [21]. At the same time, it was already shown in PIC simulations [22,23] that T e drops inside the sheath due to the deceleration caused by electric repulsion. According to the mentioned assumptions, the Vlasov equation for impinging ions reads: While energy and flux conservation can be expressed with the following equations Explicitly solving the equation (1) would require numerical techniques based for instance on PIC methods, and would lead to expensive simulations. Nonetheless, considerations about the energy distribution of particles at the wall are still possible. Similar to how the Boltzmann factor is usually derived for electron density in a retarding electric field, it can be shown that by neglecting the ⃗ v × ⃗ B term in equation (1) and by combining it with the flux and energy conservation equations (2) and (3), it follows that: where |⃗ v| is the velocity module at the SE, | ⃗ v ′ | is the velocity module at the PFM surface, and E ′ i = E i + ∆E is the total energy at the wall. Neglecting the magnetic field entails a wrong evaluation of the impinging ions' angle of incidence α (see figure 1). Indeed, as shown in [24], incident ions are still magnetized inside the sheath. Therefore, it was assumed a uniform incidence angle α = 40 • for all impinging particles, similar to what was reported in [18]. Finally, the sheath electric potential decay ϕ was approximated using the fit proposed in [25]. The data-set showed in section 3 was calculated assuming the sheath width equal to two Larmor radii of Deuterium at sonic speed multiplied by the sine of the magnetic field complementary angle (2ρ D sin(π/2 − α B )). This width well represents the potential decay described by Brooks [26]. Since for grazing magnetic fields the potential decay might be better described with the model proposed in [27], as done in [7], a comparison between the two potential decay model is discussed in section 4.4. Finally, for sputtered W atoms only electron-ion ionizations are considered, collision rates were taken from ADAS [28], and no metastable states were considered.

Boundary condition-distributions at the SE
In the present work, three different probability density functions f i (⃗ v) were tested as boundary conditions at the SE to assess the relevance of kinetic distributions in the matter of W erosion and redeposition.
1. Delta (δ) distribution. It is extensively used in the state of the art as mentioned in the introduction, it represents the average energy of impinging ions (fluid description) whereĒ SE is the average energy of ions at the SE taken from [20]. For low ionization statesĒ SE ≈ 2T i that is, the one-way heat flux associated to a zero-centered Maxwellian velocity distribution. 2. Gamma (Γ) distribution. Parallel velocity distribution described with a Gamma function while a zero centered Maxwellian is used for the perpendicular component where k 0 , θ 0 are the shape and scale parameters. The average v ∥ of impinging ions at the sheath-edge was expressed in terms of the function f(A, Z), taken from [20]. The variance was assumed to be equal to the particles' thermal velocity σ = v th = eTi mi . Where T i is the main ion temperature at the SE and m i is the incident particle mass. A gamma function was adopted because it was in good agreement with a set of incident distribution functions calculated using the PIC model shown in [29]. Essentially it was seen that the gamma function approximates the distribution function at the SE of a stream of particles originating upstream in a collisionless SOL. 3. Truncated Maxwellian (TM) distribution. Parallel velocity distribution described with a flowing Maxwellian truncated at zero for v ∥ to respect the Bohm criterion. A zero centered Maxwellian is used for the perpendicular component This case aims to model the distribution of ions accelerated in a collisional SOL. H is the Heaviside function, necessary to remove backward ions thus respecting the generalized Bohm criterion [30]. It is known that because of the so-called Knudsen layer [31] the distribution at the SE is characterized by two kinetic ionic temperatures T i ∥ and T i ⊥ . Nonetheless, in this work equipartition was imposed having Incident kinetic distributions were compared at equal mean speed and variance.

W erosion and redeposition
For a mono energetic beam of impinging ions, the energy and angular distribution of sputtered W atoms as a first approximation can be represented with a Thompson distribution as done in [20]. The angular distribution is considered cosinusoidal with respect to the polar angle (θ) and constant with respect to the azimuthal one (ϕ). The sputtered W energy-angular distribution is here shown where E ′ i is the incident projectile energy; E max is the maximum energy of sputtered W, whose formula was taken from a previous work [20]; E b = 8.67 eV is the W binding energy and dΩ = sin(θ)dθdϕ. The Thompson distribution coincides with the sputtered neutral W distribution caused by a mono energetic stream of incident particles. Normally, impinging ions' energy is taken to be equal to their average energy. Hence, their distribution is approximated with a delta function, see section 2.1. The key test assumption here is to understand if the distribution of sputtered W changes with respect to the Thompson one when the whole impinging ions energy distribution is included. For doing so it was assumed that every impinging ion with energy E ′ i greater than the sputtering threshold energy, sputters W according to a Thompson distribution Γ W +0 Th (E ′ i , E, Ω) as in equation (6). The neutral W sputtered distribution was then calculated as the convolution between the normalized Thompson distribution Γ Th,n W +0 and the impinging flux distribution, all weighted with the sputtering yield (y), as shown in equation (7) Γ is the incident flux energy distribution along z and Γ Th,n W +0 is the Thompson distribution, equation (6), normalized. Finally, α = arctan( is the impinging ions' incidence angle. As previously mentioned, the incidence angle was assumed to be fixed and equal toᾱ = 40 • .

Neutral W transport.
Once sputtered, W neutral atoms transport is described by the following Vlasov equation Neutral tungsten is subjected only to a sink term C representing ionization by electronic collisions. Therefore, C = −f W +0 n e S(n e , T e ) where n e is the local electron density and S(n e , T e ) is the ionization rate coefficient taken from ADAS [28]. As previously mentioned metastable states were neglected, thereby the ionization rate coefficient S dependance on electron density was not included. It was again assumed a uniform target and it was considered only a plasma volume in which thermal gradients are negligible. Equation (9) allows assessing the density profile of neutral W inside the sheath As expected, f W +0 (z,⃗ v) decay is related to the ionization mean free path of neutral W along the z axis: The same decay applies to the flux distribution as neutrals propagate with conserved momentum The sink term C defined in equation (9) is equal to the divergence of the neutral W flux expressed in equation (12), and hence, it is the volumic source of firstly ionized sputtered tungsten (W +1 ).

Monte Carlo simulation.
Instead of explicitly solving the Vlasov equation for W ionized states, it was decided to sample the W +1 source and push the ionized particles via Monte Carlo techniques. Firstly ionized sputtered tungsten atoms were followed within a Monte Carlo simulation to assess the probability density function of non-redeposition P nr (z, E, Ω). Assuming that the first ionization (W +0→+1 ) happens at a distance z, P nr (z, E, Ω)dEdΩ is the probability a W +1 atom with energy in (E, E + dE) and direction in dΩ has to not redeposit. The simulation time was set to five times the W +1 Larmor period as a consequence of a parametric study, see figure 2 and section 3 for more details. Particles were pushed into the monodimensional plasma background via the Boris algorithm [32] simulating the Lorentz Force. The simulation domain was limited only along the z-axis while particles were free to fly into the lateral dimensions as much as they could during the simulation period because no boundary was applied. Further ionizations in the sheath were assumed to follow an exponential distribution time as suggested in [33]. Finally, particles striking back the PFM or having velocity such that ⃗ v ·v ∥ > 0 at the end of the simulation were counted as redeposited otherwise they were considered as not redeposited. Since ionizations are an intrinsically stochastic phenomenon, more Monte Carlo simulations were performed to get a better assessment of P nr . The iteration was assumed to reach convergence with respect to the estimation of Γ net when an imposed criterion was satisfied, see equation (14) in section 3. In past works such as in [7][8][9] Monte Carlo simulations were performed to estimate the f non−redep i.e. the average probability of non-redeposition. In those works, it was defined f non−redep ≡ N non−redep /N where N is the number of simulated particles while N non−redep is the number of nonredeposited particles. In a nutshell, P nr was approximated with a binary function of zeros and ones, zero for redeposits, and one for non-redeposits. As pointed out in [9] ionizations influence f non−redep thus the net flux Γ net . For this reason, another path was followed in this paper. The probability P nr is not approximated by a binary function, the Monte Carlo simulation main output is not N non−redep but is P nr . Having the non-redeposition probability, it is possible to directly calculate the non-redeposited flux of ionized particles Γ nr W +Z as the product between P nr and the neutral W sink term C, defined in equation (9) Γ nr Adding the flux of neutral W reaching the simulation domain upper limit, considered as not redeposited, gives the total net flux Γ net . In the parameter domain studied the neutral W flux Γ nr W +0 accounted at maximum for 0.1 % of the total net flux Finally, f non−redep is defined as the ratio between the net and gross fluxes f nr ≡ Γ net /Γ gross . In a nutshell, the model here proposed differs with respect to previous works as [7][8][9] by (i) considering the impact energy distribution of impinging ions; (ii) solving analytically the W +0 decay and (iii) calculating directly Γ net instead of N non−redep .

Numerical setup
A set of Monte Carlo simulation was performed to estimate P nr for 336 different Deuterium plasma backgrounds, see parameters listed in table 1. In order to have an acceptable estimate of the probability distribution P nr (z, E, Ω)dEdΩ, and hence of Γ net , each volume dEdΩ must be sampled enough for a considered set of z. For a fixed N number of simulated particles, only η of them will be sampled in dEdΩ. Then, among the η particles sampled, a smaller number η non−redep will not be redeposited. For this reason, it was necessary to iterate the simulation n-times until the imposed convergence criterion, see equation (14), was satisfied for each simulated plasma background. Finally, P nr (z, E, Ω)dEdΩ was approximated as the corresponding fraction η/η non−redep , averaged among the whole iteration process. The simulation stopped at the nth iteration when a relative 5% difference ε = 0.05 between the net  (14) and (15), respectively A parametric scan over different combinations of total simulated time T simu and time step dt, see figure 2 was performed for a specific plasma background. As expected, the more the time step is reduced, the better the final estimate of f non−redep since the trajectories and energy balances are more accurately represented [20]. At the same time, a reduction of T simu is followed by an overestimation of f non−redep since the simulation time does not allow particles to redeposit. Finally, a time step equal to one-thousandth of the W +1 Larmor rotation period was imposed: dt = 10 −3 × τ L,W +1 and a total simulation time T simu = 5 × τ L,W +1 . The whole code was written in the programming language MATLAB ver. R2021-a. Using one core, the simulation time was around 10 min per iteration with N ∼ 2.5 × 10 5 followed particles. It is worth noting that P nr only depends on the plasma background and not on the impacting species, nonetheless, the convergence criteria does depend on it. Monte Carlo simulations were performed also to assess the IEAD discussed in section 4.2.1. In this case,  another procedure was adopted, plasma ions coming from the SOL were followed to the wall, and no P nr was estimated. It was imposed a time step dt = 10 −4 × τ L,W +1 in order to better describe the ions' angle of incidence.
An f non−redep dataset was obtained by simulating the physical sputtering of a surface made of smooth bulk W subjected to erosion by light impurities such as boron, carbon, nitrogen, oxygen (B, C, N, O), main plasma ions: D, and self-sputtering (W). Except for Deuterium, the other species were simulated up to the third ionization state. The temperature, density, and magnetic field conditions are listed in table 1. For each species, all three boundary conditions listed in section 2.1 were applied.

Model validation
The present model was validated by comparing its results with previous works. As shown in [9] if the sputtered W distribution is a Thompson distribution (see equation (6)) with E max = 2T e , f non−redep ≈ 1/2 √ χ sheath for χ sheath < 0.1, where χ sheath is the fraction of neutral W ionized outside the sheath region. It was possible to reproduce a similar cut-off energy E max = 2T e by simulating C +3 as impacting species, and using a delta distribution at the SE. As shown in figure 3 our results do not follow the f non−redep ≈ 1/2 √ χ sheath trend. This could be related to the difference in the applied potential decays and in the method used to infer f non−redep . Nonetheless, the data still correlates as a power of χ sheath times a constant. Moreover, f non−redep increases with the charge of the impinging ion. The explanation is the following: higher charge states get more energy from the plasma sheath and so they transfer more energy to the target. In turn, sputtered W has higher energy and ionization mean free path, χ sheath increases and so f non−redep . In addition to that, the model was also validated by comparing it with another previous attempt to numerically estimate W erosion, see figure 4. In [23] PIC data was used as input to the code ERO [10] to resolve the plasma sheath. Moreover, the neutral sputtered W energy and angular distribution were described again with the Thompson formula (6). Four different cases summarized in figure 5 corresponding to four different electron densities and temperatures at the sheath entrance were examined. For each case, the fraction of promptly redeposited W (f p ), was assessed four times while varying the Thompson distribution cut-off energy (E max ), see equation (6). For those cases where T e = 20 eV, the cut-off energy was set equal to 50, 100, 200, and 300 eV. Similarly, for T e = 5 eV, the cut-off energy was set equal to 10, 20, 40, and 80 eV. Again, the delta distribution was applied, and the results show a good qualitative agreement. In cases 1 and 3 the average difference is around 8 % and 5 %, figure 4. On the other hand, for T e = 5 eV (cases 2 and 4) the relative difference increases to 30 % and 100 % respectively. It is interesting to notice that except made for case 3, the present model always overestimates f p , this is probably caused by neglecting the electronic temperature drop in the sheath, thus in turn, underestimating the average mean free path of neutral W. In any case, it can be seen that neglecting the T e drop has a smaller effect if the temperature at the SE is at 20 eV. This is due to the dependence of the neutral W ionization rate coefficient S(n e , T e ) on the electron temperature. The rate coefficient from ADAS [28], is less and less dependent from T e as T e increases if T e ∈ [5, 100] eV, that essentially means that ∂ Te S(n e , T e ) drops for hotter plasma sheaths. Moreover, for low screening plasma conditions, that is, high ionization mean free path of neutral sputtered W, the fraction of sputtered particles influenced by the potential temperature drop decreases.

Neutral W non redeposition dependence on plasma parameters
As seen in section 3.1, the f non−redep value strongly depends on the fraction of neutral W ionized in the sheath. Equivalently, f non−redep is expected to depend on the average ionization mean free path of neutral W. Always considering [9] as a reference, f non−redep is correlated with the mean normalized ionization length of neutral Wλ which is defined as: where λ ion,SE is the neutral W ionization mean free path defined previously in equation (11), calculated imposing n e equal to its value at the SE. The neutral W sputtered flux at the wall Γ 0 W +0 is defined in equation (7). The normalization factor is the sheath width defined as λ sheath = 2ρ D sin(π/2 − α B ) (as in [25]). The parameterλ was fitted with a scaling law to relate it to plasma parameters, see equations (17) and (18). The fit here proposed has an accuracy within ±40 percent. It is valid for magnetic field angles α B ∈ [2, 85] deg, and for impinging ion mass ranging from D to W, see section 3 to see all the simulated conditions. Equation (17) where the units of measure are: 10 19 ×m −3 , eV, T, deg, amu, 1. It is worth noting that the ionization mean free path depends roughly on the inverse of plasma pressure (n e T e ) and that the impinging ion mass dependence is very weak, as in [34]. f non−redep trends as function ofλ were compared using both the fluid (δ) and kinetic (Γ,TM) boundary conditions. From figures 6(a) and (b) the gap between Γ and TM is small, suggesting that a truncated Maxwellian (TM) in the parallel direction is already a good approximation for impinging ions distribution at the SE. Actually, this is not always true, since for light plasma ions such as Deuterium the gap between Γ and TM can be large, see figure 10. Due to its low mass, Deuterium distributions have a high variance. This makes Deuterium more sensitive to the chosen boundary condition. Moreover, it was seen that choosing a Gamma distribution always results in a higher f non−redep due to its longer high-energy tail. Finally, adopting the whole kinetic energy distribution is important only for highly screened plasmasλ ≪ 1, and for impinging ions with low ionization states. For N +1 using a delta function leads to an underestimation of f non−redep up to 10 times whenλ is at its minimum value. On the other hand, for N +3 the maximum difference is around 40 percent.

Magnetic field configuration influence on W redeposition.
Variations in the magnetic field intensity and angle of incidence affect the redeposited fraction. For instance, an increase in the intensity results in a decrease of the Larmor radius ρ D thus in a decrease of the plasma sheath thickness or equivalently, an increase of theλ factor, see [9] and equation (17). Similarly, as the angle of incidence α B increases, the sheath thickness decreases. In both cases, by reducing the sheath width there is an increase inλ, hence an increase of the net W source. The increase in f non−redep as | ⃗ B| or α B increases depends onλ. In general, the lower the electron density n e and temperature T e (higherλ) the less f non−redep is sensitive to plasma conditions, see figures 7(a) and (b). Actually, forλ ≫ 1, since most particles ionize outside the sheath, also W falls in the redeposition ballistic regime described by Fussman et al [6]. The relationship betweenλ and f non−redep is not always perfectly monotonically increasing in our database. As a consequence, especially when it is calculated the ratio between two small f non−redep values some noise may appear as in figure 7(b) where an outlier is found. Moreover, what is shown in the figure 7(b), demonstrates that in order to compare W erosion studies done in different machines it is necessary to scale the data according to the difference in intensity between the applied magnetic fields.

Effective sputtering yield
The sputtering yield y(E ′ i , α) is a function of the energy and angle of impact of the incident particles. It quantifies the number of eroded particles per incident particle and is fundamental for the quantification of the gross eroded flux Γ gross = Γ inc ×Ȳ. The sputtering yield formula for the database is that proposed by Garcia in [35]. As shown in [36], and done in [14] the sputtering yield y shall be integrated with the distribution function of impinging particles at the wall f i (E ′ , α ′ ) resulting in an effective sputtering yieldȲȲ The assumption of incident distribution equal to a Truncated Maxwellian or a gamma function thus results in a different estimate ofȲ than in the fluid case where a δ distribution is instead considered. Fixing the impinging ion angle of incidence α = 40 • remove the angular dependence, thus,Ȳ depends only on the impact energy. The total impact energy in turn mainly depends on the electronic temperature T e , the ratio τ = T i /T e in the sheath, and the impinging ion ionization state. It was possible to expressȲ as a second-order polynomial function of the sputtered W temperature, which was defined as Where, for the light impurities B, C, N, O, the temperature T W +0 can be approximated by a scaling law, T w ≈ ξ 1 T ξ2 e Z ξ3 , with fitting errors less than 1 percent if τ = 1. The C and ξ coefficients of equation (20) slightly depend on the ions mass (m i ). As shown in figure 8, considering the incident energy distribution can bring up to a difference between 4 and 5 times in the sputtering yield of firstly ionized light impurities. On the other hand, for the second and third ionization states the difference is around 10% ÷ 30%. The energy distribution seems to have in general a negligible impact on more than once ionized particles. The ratio between the effective sputtering yield calculated with a delta and a TM functions according to equation (19) decreases as the electronic temperature and the impinging ion charge increases. Essentially, the higher the energy gain inside the sheath the smaller the difference. As opposed to f non−redep , an increase in T e reduces the difference between the δ and TM distributions. This shows that for W redeposition (f non−redep ) the increase in ionization rate is faster than the increase in sputtered W ionization mean free path for increasing T e .

IEAD, the impact energy and angular distribution.
In section 4.1 was shown the difference in f non−redep between fluid (δ) and kinetic (Γ or TM) sputtering description while considering a fixed impact angle. In the following section, we temporarily relax this assumption by adopting a Monte Carlo simulation to estimate f i (E ′ i , α) hence the impinging ions energy and angular distribution (or IEAD). A TM distribution of impinging ions (see section 2.1) was sampled at the sheath entrance. Particles were pushed using the Boris [32] algorithm till the PFM surface. As a preliminary study it was simulated a stream of O +1 immersed in high screening plasma conditions (λ ≪ 1). First, we noted (see figure 9(b)) that the description of the energy distribution given by the equations (1)-(3) matches the Monte Carlo simulation estimate for O +1 and T e = 15 eV. Moreover, as shown in figure 9(a), another factor of six rose up between f non−redep estimates with a fixed angle of incidence or by considering the whole IEAD whenλ ≪ 1. Finally, a 30% to 60% increase in the sputtering yield was calculated. Again, this was verified only for O +1 , further studies should be done. However, a decrease in the f non−redep ratio shown in figure 9(a) is envisaged for higher charge states. Moreover, because of roughness [37], the sputtering yield angular dependence decreases reducing the effect related to the angular distribution of impinging ions.

Plasma ions erosion-Deuterium
By considering the whole energy distribution of impacting ions, also those species whose average energy was below the sputtering threshold can still erode because of their high-energy tail. It was possible to calculate the f non−redep caused by D +1 in any conditions, see figure 10. Since D +1 is one of the most abundant species inside tokamaks operating with it as a fuel, assessing its relative contribution to W erosion could be of interest. First of all, if the erosion provoked by different impinging impurities is considered a linear process, the total flux Γ net is given by the sum of all contributing species (21) where Γ inc is the total incident flux, Cc is the impurity concentration, andȲ net ≡ f non−redepȲ is the net sputtering yield. If we consider every impurity having the same concentration, their contribution can be compared looking at the net sputtering yield. In order to make an accurate comparison, it was decided to use the experimental sputtering yield values taken from Behrisch and Eckstein ( [38]). Given the lack of data, the comparison was made only between Deuterium, Nitrogen, and Tungsten. As mentioned in section 4.1, the net sputtering yield of Deuterium strongly depends on the boundary condition. From figure 10, it is clear that (Ȳ net D +1 ) is always the smallest. Nonetheless, the net flux also depends on the impurity concentration. If impurities concentration with ionization state Z = 1 is around ∼ 10 −2 , Deuterium and impurities'Ȳ net is comparable. For instance, as shown in figure 10, Deuterium Y net is expected to be less than 100 times smaller whenλ > 0.3. Hence, there is a set of plasma conditions for which Deuterium contribution to the total net flux is not negligible.

Potential drop model uncertainty
Once plasma conditions ( ⃗ B, n e , T e , T i ) and the incident ion mass and charge (m i , Z i ) are fixed, the fraction f non−redep depends mainly on the sheath thickness. The potential drop is usually described with the Brooks or Stangeby [27] models. The difference between the two models concerns the sheath thickness. In the electric potential drop fit proposed in [25] and used here, to reproduce Stangeby's electric potential profiles one has to enlarge the sheath thickness from 2ρ D sin(π/2 − α B ) for Brooks model to 3ρ D sin(π/2 − α B ) for Stangeby's. As shown in [7], the uncertainty on the model leads to a difference in the estimation of f non−redep asλ falls. Consequently, the difference in f non−redep made with a fluid or kinetic model also increases. Using the potential drop proposed by Stangeby, the f non−redep related to incident light impurities ionized once with a fluid model can be three orders of magnitude smaller than the one calculated considering the incident energy distribution, see figure 11(a). Nonetheless, the kinetic correction saturates quickly with respect to the impinging ion charge. As shown in figure 11(b), the kinetic correction becomes less important resulting in a factor of two difference whenλ ≪ 1. Finally, there is also an uncertainty related to the ionization rates. Evidence have been gathered as in [39] according to which neutral W metastable states should be considered to better estimate W ionization coefficients. This would in theory reduce the average ionization length of neutral tungsten, thusλ.

Discussion and conclusions
In this work, the fraction of non-redeposited W was calculated for a wide range of simplified 1D plasma conditions by a semi-analytical model based on Monte Carlo techniques. The model was benchmarked with previous attempts to numerically assess net erosion and it showed a good agreement in the parametrical range of this study. The results respected the expected physical trends. It was seen, for example, that the fraction of non-redeposited W depends mainly on the competition between the average ionization length of the sputtered neutral W and the sheath thickness. The scan of 336 plasma backgrounds and the sputtering induced by different impurities ranging from Deuterium to Tungsten showed that the kinetic description of sputtering accounting for the incident energy distribution, with fixed impact angle, has a nonnegligible effect only for impurities ionized once and for sputtering produced by plasma ions such as Deuterium (sections 4.1 and 4.3). The explanation is the following: the refined distribution of sputtered neutral W following the kinetic description is nothing more than a weighted average between the incident distribution, the sputtering yield, and the Thompson distribution. It is thus intuited that if the sputtering yield is constant, the distribution of the sputtered neutral W tends to the Thompson one caused by a mono energetic beam of particles whose energy is equal to the distribution average energy. The sputtering yield variation decreases rapidly when the energy of the impinging ion is greater than its minimum sputtering energy; E ′ i > E min . In general, the sputtering yield tends to a constant value for E i + ∆E >> E min . Since ∆E ≈ ∆ϕ Z the impact energy of the incident ions is directly proportional to their ionization state Z and to the sheath electron temperature T e . According to the model discussed here, when T e ⩾ 15eV, already for Z = 2 the sputtering yield is roughly constant. Thus, particles in the fast tails of the distribution will contribute about the same as the less energetic particles because of an overkill effect where the excess energy is not increasing much the sputtering yield. However, even higher ionization states are expected to benefit from kinetic correction if the electronic temperature is very low. It is possible to demonstrate the above by taking the probability density function (PDF) associated with the equation (7) Γ 0,n W +0 (E, Ω) =´E If the incident particle energy distribution is such that the sputtering yield is approximately constant and thus independent of the incident energy, it can be simplified. Therefore, the PDF tends to that produced by a delta function centered in the mean value, i.e. to the fluid value Γ 0,n W +0 (E, Ω) =´E WithĒ ′ i being the average incident energy. So, for a fixed angle of incidence, the sputtering yield converges to a constant for most species. Despite that, the convergence actually depends on the angle of incidence as well. For example, for some angles, the convergence is reached for higher energies and the kinetic correction might be more pronounced, (section 4.2.1). In this work, the kinetic correction regards only the neutral sputtered W high-energy tail and not its angular distribution. The importance of this correction depends on how much the high-energy tail population increases by it and on plasma screeningλ. If plasma screening is low (λ ≳ 1) most of the particles are not redeposited and the population fraction populating the tails is negligible if compared to the whole nonredepositing population. At the same time, the tail population is less affected by the kinetic correction for ionization states greater than one, see sections 4.1 and 4.4. On top of that, it shall be kept in mind that by considering reflections the transferred energy to the target decreases. Therefore, the sputtering yield is more sensitive to the energy distribution of impinging ions thus the kinetic correction could be more important. The sputtering yield energy derivative is indeed higher for lower impinging energies. On the other hand, the electron temperature drop inside the sheath tends to diminish the kinetic correction increasingλ. Lastly, the material roughness could also vary the sputtering yield sensitivity to the incident energy. For instance, it was seen that the sputtering yield angular dependence flattens out if the roughness is included [37]. The overall balance between these effects mentioned should be analyzed in future proceedings. Finally, in this work, it was possible to compare erosion from Deuterium with that produced by light and heavy impurities for T e ∈ [15,30] eV. It was seen that the contribution of Deuterium is not always negligible if compared to once-ionized impurities. In the future, it is planned to use the data produced by the model described in this paper to generate a fitting function capable to predict net W sources, that is the effective sputtering yield. This could in turn be used in fluid plasma transport codes to auto-consistently estimate W sources.

Data availability statement
The data cannot be made publicly available upon publication because they are not available in a format that is sufficiently accessible or reusable by other researchers. The data that support the findings of this study are available upon reasonable request from the authors.