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

L. Cappelli, N. Fedorczak, J. P. Gunn, S. Di Genova, J. Guterl, E. Serre

► To cite this version:

L. Cappelli, N. Fedorczak, J. P. Gunn, S. Di Genova, J. Guterl, et al. Study of the erosion and redeposition of W considering the kinetic energy distribution of incident ions through a semi-analytical model. Plasma Physics and Controlled Fusion, 2023, 65 (9), pp.095001. 10.1088/1361-6587/ace282. hal-04190861

HAL Id: hal-04190861

https://hal.science/hal-04190861v1

Submitted on 30 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Study of the erosion and redeposition of W
considering the kinetic energy distribution of
incident ions through a semi-analytical model

L Cappelli<sup>1,2,*</sup>, N Fedorczak<sup>2</sup>, J P Gunn<sup>2</sup>, S Di Genova<sup>1</sup>, J Guterl<sup>3</sup>, E Serre<sup>1</sup>

<sup>1</sup> M2P2 Aix-Marseille Univ, CNRS, Centrale Marseille, 13013 Marseille, France

<sup>2</sup> CEA, IRFM, F-13108 Saint-Paul-lez-Durance, France

<sup>3</sup> General Atomics, PO Box 85608, San Diego, CA 92186-5608, USA

Abstract

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.

**Keywords:** net sources, plasma-wall interactions, tungsten, erosion, redeposition, kinetic dis-tributions, Monte Carlo.

# 1 Introduction

In the future, magnetic nuclear fusion reactors might rely on metallic plasma-facing components (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

<footnote>*E-mail: luca.cappelli@univ-amu.fr</footnote>
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’ (Γ<sub>gross</sub>), which can be expressed as the
incident flow (Γ<sub>inc</sub>) times an effective sputtering yield (Y). 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<sub>p</sub> (non-prompt fraction f<sub>np</sub> = 1 − f<sub>p</sub>), and it was
used to adjust the sputtering yield redefined as ‘net’ sputtering yield. This led to a distinction
between the gross flux Γ<sub>gross</sub> = Γ<sub>inc</sub> × Y and the net flux Γ<sub>net</sub> = Γ<sub>gross</sub> × f<sub>np</sub>. It was shown by
[7] that the definition of f<sub>np</sub> adopted in [6] was not appropriate for high-Z PFMs. In particular,
through a set of local Monte Carlo simulations, they noted that f<sub>np</sub> 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<sub>np</sub> was redefined as f<sub>non-redep</sub>, 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<sub>non-redep</sub>, 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<sub>i</sub> is
expressed as a function of integrated or ‘fluid’ plasma quantities, e.g. E<sub>i</sub> = 2T<sub>i</sub> + Δφ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<sub>non-redep</sub> 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.

## 2 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(\vec{v})$ at the sheath entrance (SE) accelerated inside the plasma sheath. Ions are immersed in a uniform magnetic field B oblique with angle $\alpha_B$ respect to an ideal flat surface (no $\nabla \vec{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 \approx 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:

$$v_z \partial_z f_i + \frac{Ze}{m} (\vec{E} + \vec{v} \times \vec{B}) \cdot \nabla_v f_i = 0 \quad (1)$$

While energy and flux conservation can be expressed with the following equations.

$$f_i(\vec{v}) v_z d\vec{v} = f_i(\vec{v}') v_z' d\vec{v}' \quad (2)$$

$$v_z'^2 = v_z^2 + \frac{2Ze\Delta\phi}{m} \quad (3)$$

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 $\vec{v} \times \vec{B}$ term in equation (1) and by combining it with the flux and energy conservation equations (2), (3), it follows that:

$$f_i(|\vec{v}|) = f_i(|\vec{v}'|) \quad (4)$$

$$f_i(|\vec{v}'|) d|\vec{v}'| = f_i(E'_i) dE' \quad (5)$$

Where $|\vec{v}|$ is the velocity module at the SE, $|\vec{v}'|$ is the velocity module at the PFM surface, and $E'_i = E_i + \Delta E$ is the total energy at the wall. Neglecting the magnetic field entails a wrong evaluation of the impinging ions' angle of incidence $\alpha$ (see figure 1). Indeed, as shown
<figure><img src="image_2.png" /><figcaption>Figure 1: Coordinate systems, top left for incident ions, bottom right for sputtered W. Particles and magnetic field $\vec{B}$ etching with angles $\alpha$ and $\alpha_B$ respectively.</figcaption></figure>

in [24], incident ions are still magnetized inside the sheath. Therefore, it was assumed a uniform incidence angle $\alpha = 40^{\circ}$ for all impinging particles, similar to what was reported in [18]. Finally, the sheath electric potential decay $\phi$ 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 times the multiplied by the sine of the magnetic field incidence angle $(2\rho_D\sin(\pi/2 - \alpha_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.

## 2.1 Boundary condition - distributions at the Sheath Entrance

In the present work, three different probability density functions $f_i(\vec{v})$ were tested as boundary conditions at the Sheath Entrance (SE) to assess the relevance of kinetic distributions in the matter of W erosion and redeposition.

1. Delta ($\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).

$$f_i(E_i) = \delta(E_i - \bar{E}_{SE})$$

Where $\bar{E}_{SE}$ is the average energy of ions at the SE taken from [20]. For low ionization states $\bar{E}_{SE} \approx 2T_i$ that is, the one-way heat flux associated to a zero-centered Maxwellian velocity distribution.

2. Gamma ($\Gamma$) distribution. Parallel velocity distribution described with a Gamma function while a zero centered Maxwellian is used for the perpendicular component.

$$f_i(v_{\parallel}, v_{\perp}) = \frac{v_{\parallel}^{k_0-1}}{\Gamma(k_0)\theta^{k_0}} e^{-\frac{v_{\parallel}}{\theta_0}} \sqrt{\frac{B_{\perp}}{\pi}} e^{-B_{\perp}v_{\perp}^2}$$

where $k_0, \theta_0$ are the shape and scale parameters. The average $v_{\parallel}$ 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 $\sigma = v_{th} = \sqrt{\frac{eT_i}{m_i}}$. 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_\parallel$ to respect the Bohm criterion. A zero centered Maxwellian is used for the perpendicular component.

$$f_i(v_\parallel, v_\perp) = \sqrt{\frac{B_\parallel}{\pi}} e^{-B_\parallel v_\parallel^2} \sqrt{\frac{B_\perp}{\pi}} e^{-B_\perp v_\perp^2} H(v_\parallel)$$

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_\parallel^i$ and $T_\perp^i$. Nonetheless, in this work equipartition was imposed having $T_\parallel^i = T_\perp^i$.

Incident kinetic distributions were compared at equal mean speed and variance.

## 2.2 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 ($\theta$) and constant with respect to the azimuthal one ($\phi$). The sputtered W energy-angular distribution is here shown.

$$\Gamma_{W+0}^{Th}(E'_i, E, \Omega)dEd\Omega \propto \frac{1}{2\pi(E+E_b)^3} \left[ 1 - \sqrt{\frac{E+E_b}{E_{max}(E_i)+E_b}} \right] \cos(\theta)dEd\Omega \quad (6)$$

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\Omega = \sin(\theta)d\theta d\phi$. 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 $\Gamma_{Th}^{W+0}(E'_i, E, \Omega)$ as in equation (6). The neutral W sputtered distribution was then calculated as the convolution between the normalized Thompson distribution $\Gamma_{W+0}^{Th,n}$ and the
impinging flux distribution, all weighted with the sputtering yield (y), as shown in equation (7).

$$ \Gamma_{W+0}^{0}(E, \Omega) = \int_{E'_{i}} \Gamma_{i}(E'_{i})y(E'_{i}, \bar{\alpha})\Gamma_{W+0}^{Th,n}(E'_{i}, E, \Omega)dE'_{i}dEd\Omega \quad (7) $$

$$ \Gamma_i(\vec{v}) = f_i(\vec{v})v_z \quad (8) $$

$\Gamma_i(E'_i)$ is the incident flux energy distribution along z and $\Gamma_{W+0}^{Th,n}$ is the Thompson distribution, equation (6), normalized. Finally, $\alpha = \arctan\left(\frac{\sqrt{v_x^2+v_y^2}}{v_z}\right)$ is the impinging ions' incidence angle. As previously mentioned, the incidence angle was assumed to be fixed and equal to $\bar{\alpha} = 40^{\circ}$.

### 2.2.1 Neutral W transport

Once sputtered, W neutral atoms transport is described by the following Vlasov equation.

$$ v_z \partial_z f_{W+0} = C \quad (9) $$

Neutral tungsten is subjected only to a sink term C representing ionization by electronic collisions. Therefore, $C = -f_{W+0}n_eS(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 dependence 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 $f_{W+0}(z, \vec{v})$, see equation (10).

$$ f_{W+0}(z, \vec{v}) = f_{W+0}(0, \vec{v})e^{-\int_{0}^{z} \frac{n_e(z)S(T_e)}{v_z} dz} \quad (10) $$

As expected, $f_{W+0}(z, \vec{v})$ decay is related to the ionization mean free path of neutral W along the z axis:

$$ \lambda_{ion,z} = \frac{v_z}{n_e(z)S(T_e)} \quad (11) $$

The same decay applies to the flux distribution as neutrals propagate with conserved momentum.

$$ \Gamma_{W+0}(z, \vec{v}) = \Gamma_{W+0}^{0}(\vec{v})e^{-\int_{0}^{z} \frac{dz}{\lambda_{ion,z}}} \quad (12) $$

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}$).
### 2.2.2 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, \Omega)$. Assuming that the first ionization ($W^{+0} \to ^{+1}$) happens at a distance $z$, $P_{nr}(z, E, \Omega)dEd\Omega$ is the probability a $W^{+1}$ atom with energy in ($E, E + dE$) and direction in $d\Omega$ has to not redeposit. The simulation time was set to 5 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 $\vec{v} \cdot \hat{v}_{\parallel} > 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 $\Gamma_{net}$ when an imposed criterion was satisfied, see equation (14) in section 3. In past works such as in [7], [9] and [8] 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} \equiv N_{non-redep}/N$ where $N$ is the number of simulated particles while $N_{non-redep}$ is the number of non-redeposited 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 $\Gamma_{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 $\Gamma_{W+z}^{nr}$ as the product between $P_{nr}$ and the neutral W sink term $C$, defined in equation (9).

$$ \Gamma_{W+z}^{nr} = \int C(z, E, \Omega)P_{nr}(z, E, \Omega)dzdEd\Omega \quad (13) $$

Adding the flux of neutral W reaching the simulation domain upper limit, considered as not redeposited, gives the total net flux $\Gamma_{net}$. In the parameter domain studied the neutral W flux $\Gamma_{W+0}^{nr}$ accounted at maximum for 0.1 % of the total net flux.

$$ \Gamma_{net} = \Gamma_{W+z}^{nr} + \Gamma_{W+0}^{nr} $$

Finally, $f_{non-redep}$ is defined as the ratio between the net and gross fluxes $f_{nr} \equiv \Gamma_{net}/\Gamma_{gross}$. In a nutshell, the model here proposed differs with respect to previous works as [9], [7] and [8] by i) considering the impact energy distribution of impinging ions; ii) solving analytically the $W^{+0}$ decay and iii) calculating directly $\Gamma_{net}$ instead of $N_{non-redep}$.
## 3 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, \Omega)dEd\Omega$, and hence of $\Gamma_{net}$, each volume $dEd\Omega$ must be sampled enough for a considered set of $z$. For a fixed $N$ number of simulated particles, only $\eta$ of them will be sampled in $dEd\Omega$. Then, among the $\eta$ particles sampled, a smaller number $\eta_{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, \Omega)dEd\Omega$ was approximated as the corresponding fraction $\eta/\eta_{non-redep}$, averaged among the whole iteration process. The simulation stopped at the nth iteration when a relative 5 percent difference $\epsilon = 0.05$ between the net fluxes estimated at the nth and nth-1 iterations was attained. The convergence criteria and the non-redeposition probability updated to the nth iteration are here shown in equations (14) and (15), respectively.

Table 1: plasma parameters used to simulate W erosion.

<table><tr><td>$T_i = T_e$</td><td>15, 20, 25, 30</td><td>eV</td></tr><tr><td>$n_e$</td><td>5e18, 1e19, 5e19, 1e20</td><td>m<sup>-3</sup></td></tr><tr><td>$B$</td><td>2, 3.7, 5</td><td>T</td></tr><tr><td>$\alpha_B$</td><td>2, 4, 8, 20, 40, 75, 85</td><td>deg</td></tr></table>

$$ \left| \frac{\Gamma_{net}^n}{\Gamma_{net}^{n-1}} - 1 \right| \le \epsilon \qquad (14) $$

$$ P_{nr}^{n}(z, E, \Omega)dEd\Omega \approx \frac{1}{n} \sum_{i}^{n} \left( \frac{\eta_{non-redep}}{\eta} \right)_{i} \qquad (15) $$

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} \times \tau_{L,W+1}$ and a total simulation time $T_{simu} = 5 \times \tau_{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 minutes per iteration with $N \sim 2.5 \times 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} \times \tau_{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,
<figure><img src="image_3.png" /><figcaption>Figure 2: Parametric scan of different *dt*, *T*<sub>simu</sub> combinations normalized with respect to the *W*<sup>+1</sup> Larmor period. Each tested point is represented with a red square, while the chosen combination is highlighted with a star. The interpolated 2D map refers to the following case: *n*<sub>e</sub> = 1e19 m<sup>-3</sup>, *T*<sub>i</sub> = *T*<sub>e</sub> = 20 eV, *B* = 3.7 T, *α*<sub>B</sub> = 2 deg, impinging species: *O*<sup>+1</sup>, δ model.</figcaption></figure>

and magnetic field conditions are listed in table 1. For each species, all 3 boundary conditions listed in section 2.1 were applied.

## 3.1 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} \approx 1/2\sqrt{\chi_{sheath}}$ for $\chi_{sheath} < 0.1$, where $\chi_{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

<table><caption>Figure 3: List of investigated plasma $n_e$ and $T_e$ at the SE. Image adapted from [23]</caption><thead><tr><th>Case #</th><th>T<sub>e</sub> (eV) at sheath entrance</th><th>n<sub>e</sub> (cm<sup>-3</sup>) at sheath entrance</th></tr></thead><tbody><tr><td>1</td><td>20</td><td>6 &times; 10<sup>13</sup></td></tr><tr><td>2</td><td>5</td><td>6 &times; 10<sup>13</sup></td></tr><tr><td>3</td><td>20</td><td>5 &times; 10<sup>11</sup></td></tr><tr><td>4</td><td>5</td><td>5 &times; 10<sup>11</sup></td></tr></tbody></table>

As shown in figure 4 our results do not follow the $f_{non-redep} \approx 1/2\sqrt{\chi_{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 $\chi_{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, $\chi_{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 5. 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 3 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 5. 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 \in [5, 100]$ eV, that essentially means that $\partial_{T_e}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.

<figure><img src="image_4.png" /><figcaption>Figure 4: Incident $\delta$ distribution. $f_{non-redep}$ trends as function of $\chi_{sheath}$, $\alpha_B$ = 2°.</figcaption></figure>

# 4 Results

## 4.1 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
<figure><img src="image_5.png" /><figcaption>Figure 5: Numerical estimates of prompt redeposition $f_p$ for the plasma conditions from table 1 in dependence on the cut-off energy. Comparison between ERO-PIC and present work results, a) $T_e = 20$ eV, b) $T_e = 5$ eV.</figcaption></figure>

mean normalized ionization length of neutral W $\hat{\lambda}$ which is defined as:

$$ \hat{\lambda} = \frac{1}{\Gamma_{gross} \lambda_{sheath}} \int \lambda_{ion,SE} \Gamma_{W+0}^{0}(E, \Omega) dEd\Omega \quad (16) $$

Where $\lambda_{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 sheath entrance (SE). The neutral W sputtered flux at the wall $\Gamma_{W+0}^0$ is defined in equation (7). The normalization factor is the sheath width defined as $\lambda_{sheath} = 2\rho_D \sin(\pi/2 - \alpha_B)$ (as in [25]). The parameter $\hat{\lambda}$ was fitted with a scaling law to relate it to plasma parameters, see equation (17), (18). The fit here proposed has an accuracy within ± 40 percent. It is valid for magnetic field angles $\alpha_B \in [2, 85] deg$, and for impinging ion mass ranging from D to W, see section 3 to see all the simulated conditions. Equation (17) was calculated assuming a fluid plasma description.

<figure><img src="image_6.png" /><figcaption>Figure 6: Differences in $f_{non-redep}$ trends for the three different boundary conditions $\delta$, $\Gamma$, TM. Grazing magnetic field, $\alpha_B = 2^\circ$. a) $N^{+1}$, b) $N^{+3}$</figcaption></figure>
$$ \hat{\lambda} \approx g(\alpha_B) n_e^{-1} T_e^{-0.88} B m_i^{0.01} Z_i^{0.2} \quad (17) $$

$$ g(\alpha_B) = 9.3 \times 10^{-6} \alpha_B^4 - 1.2 \times 10^{-3} \alpha_B^3 + 4.7 \times 10^{-2} \alpha_B^2 - 0.56 \alpha_B + 5.4 \quad (18) $$

Where the units of measure are: $10^{19} \times 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 $\hat{\lambda}$ were compared using both the fluid ($\delta$) and kinetic (Γ,TM) boundary conditions. From figures 6a and 6b 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 $\hat{\lambda} \ll 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 $\hat{\lambda}$ is at its minimum value. On the other hand, for $N^{+3}$ the maximum difference is around 40 percent.

### 4.1.1 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 $\rho_D$ thus in a decrease of the plasma sheath thickness or equivalently, an increase of the $\hat{\lambda}$ factor, see [9] and equation (17). Similarly, as the angle of incidence $\alpha_B$ increases, the sheath thickness decreases. In both cases, by reducing the sheath width there is an increase in $\hat{\lambda}$, hence an increase of the net W source. The increase in $f_{non-redep}$ as $|\vec{B}|$ or $\alpha_B$ increases depends on $\hat{\lambda}$. In general, the lower the electron density $n_e$ and temperature $T_e$ (higher $\hat{\lambda}$) the less $f_{non-redep}$ is sensitive to plasma conditions, see figures. 7a and 7b. Actually, for $\hat{\lambda} \gg 1$, since most particles ionize outside the sheath, also W falls in the redeposition ballistic regime described by Fussman [6]. The relationship between $\hat{\lambda}$ 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 7b where an outlier is found. Moreover, what is shown in the figure 7b, 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.

## 4.2 Effective sputtering yield

The sputtering yield $y(E'_i, \alpha)$ 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
<figure><img src="image_7.png" /><figcaption>Figure 7: Differences in $f_{non-redep}$ trends for the three different boundary conditions $\delta$, $\Gamma$, TM. a) increase in $f_{non-redep}$ going from $\alpha_B = 4^\circ$ to $\alpha_B = 85^\circ$, b) increase in $f_{non-redep}$ going from $B = 2 T$ to $B = 5 T$ with fixed $\alpha_B = 4^\circ$.</figcaption></figure>

for the quantification of the gross eroded flux $\Gamma_{gross} = \Gamma_{inc} \times \bar{Y}$. The sputtering yield formula for the database is that proposed by Garcia in [35]. As shown in [36], and done in [14] $y$ shall be integrated with the distribution function of impinging particles at the wall $f_i(E', \alpha')$ resulting in an effective sputtering yield $\bar{Y}$.

$$ \bar{Y} = \frac{1}{\Gamma_{inc}} \iint \Gamma_i(E'_i, \alpha) y(E'_i, \alpha) dE'_i d\alpha \quad (19) $$

The assumption of incident distribution equal to a Truncated Maxwellian or a gamma function thus results in a different estimate of $\bar{Y}$ than in the fluid case where a $\delta$ distribution is instead considered. Fixing the impinging ion angle of incidence $\alpha = 40^\circ$ remove the angular dependence, thus, $\bar{Y}$ depends only on the impact energy. The total impact energy in turn mainly depends on the electronic temperature $T_e$, the ratio $\tau = T_i/T_e$ in the sheath, and the impinging ion ionization state. It was possible to express $\bar{Y}$ as a second-order polynomial function of the sputtered W temperature, which was defined as $T_{W+0} \equiv \Gamma_{gross}^{-1} \int \Gamma_{W+0}^0(E, \theta) EdEd\theta$.

$$ \bar{Y} \approx C_1 T_{W+0}^2 + C_2 T_{W+0} + C_3 \quad (20) $$

Where, for the light impurities $B, C, N, O$, the temperature $T_{W+0}$ can be approximated by a scaling law, $T_w \approx \xi_1 T_e^{\xi_2} Z^{\xi_3}$, with fitting errors less than 1 percent if $\tau = 1$. The $C$ and $\xi$ 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 $\div$ 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 $\delta$ 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$.

<figure><img src="image_8.png" /><figcaption>Figure 8: Ratio of $\bar{Y}$ calculated using ($\delta$) and TM as incident distributions. Results are shown for boron in its first three ionization states, $B^{+1}$, $B^{+2}$, $B^{+3}$.</figcaption></figure>

### 4.2.1 IEAD, the Impact Energy and Angular Distribution

<figure><img src="image_9.png" /><figcaption>Figure 9: a) Incident $O^{+1}$ on W. Ratio of kinetic $f_{non-redep}$ both using a TM distribution with a fixed angle of incidence $\bar{\alpha} = 40^\circ$ and the whole IEAD. b) Energy distribution at the target estimated for the two cases.</figcaption></figure>

In section 4.1 was shown the difference in $f_{non-redep}$ between fluid ($\delta$) 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, \alpha)$ 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 ($\hat{\lambda} \ll 1$). First, we noted (see figure 9b) that the description of the energy distribution given by the equations (1 to 3) matches the Monte Carlo simulation estimate for $O^{+1}$ and $T_e = 15$ eV. Moreover, as shown in figure 9a, another factor of 6 rose up between $f_{non-redep}$ estimates with a fixed angle of incidence or by considering the whole IEAD when $\hat{\lambda} \ll 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 9a 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.

## 4.3 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 $\Gamma_{net}$ is given by the sum of all contributing species.

$$ \Gamma_{net} = \Gamma_{inc} \sum_{i} (CcY^{\bar{net}})_{i} \qquad (21) $$

Where $\Gamma_{inc}$ is the total incident flux, $Cc$ is the impurity concentration, and $Y^{\bar{net}} \equiv f_{non-redep}Y$ 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 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 $(Y^{\bar{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 $\sim 10^{-2}$, Deuterium and impurities' $Y^{\bar{net}}$ is comparable. For instance, as shown in figure 10, Deuterium $Y^{\bar{net}}$ is expected to be less than 100 times smaller when $\hat{\lambda} > 0.3$. Hence, there is a set of plasma conditions for which Deuterium contribution to the total net flux is not negligible.

## 4.4 Potential drop model uncertainty

Once plasma conditions $(\vec{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
<figure><img src="image_10.png" /><figcaption>Figure 10: $\bar{Y} \times f_{non-redep}$ trends with respect to $\hat{\lambda}$ for $D^{+1}, N^{+1}, W^{+1}$. Incident TM and $\Gamma$ distributions, $\alpha_B = 2^\circ$ and $T_e = T_i = 30eV$.</figcaption></figure>

from $2\rho_D\sin(\pi/2 - \alpha_B)$ for Brooks model to $3\rho_D\sin(\pi/2 - \alpha_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 $\hat{\lambda}$ 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 3 orders of magnitude smaller than the one calculated considering the incident energy distribution, see figure 11a. Nonetheless, the kinetic correction saturates quickly with respect to the impinging ion charge. As shown in figure 11b, the kinetic correction becomes less important resulting in a factor of 2 difference when $\hat{\lambda} \ll 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 $\hat{\lambda}$.

# 5 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
<figure><img src="image_11.png" /><figcaption>Figure 11: Stangeby's potential drop; differences in $f_{non-redep}$ trends for the three different boundary conditions $\delta$, $\Gamma$, TM for a grazing magnetic field, $\alpha_B = 2^\circ$. a) $N^{+1}$, b) $N^{+2}$</figcaption></figure>

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 + \Delta E \gg E_{min}$. Since $\Delta E \approx \Delta\phi 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 \ge 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).

$$ \Gamma_{W+0}^{0,n}(E, \Omega) = \frac{\int_{E'_i} \Gamma_i(E'_i) y(E'_i, \bar{\alpha}) \Gamma_{W+0}^{Th,n}(E'_i, E, \Omega) dE'_i dE d\Omega}{\int_{E'_i} \Gamma_i(E'_i) y(E'_i, \bar{\alpha}) dE'_i} $$

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.

$$ \Gamma_{W+0}^{0,n}(E, \Omega) = \frac{\int_{E'_i} \Gamma_i(E'_i) \Gamma_{W+0}^{Th,n}(E'_i, E, \Omega) dE'_i dE d\Omega}{\int_{E'_i} \Gamma_i(E'_i) dE'_i} = \Gamma_{W+0}^{0,n}(\bar{E}'_i, E, \Omega) $$

With $\bar{E}'_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 $\hat{\lambda}$. If plasma screening is low ($\hat{\lambda} \gtrsim 1$)
most of the particles are not redeposited and the population fraction populating the tails is
negligible if compared to the whole non-redepositing 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.
On the other hand, the electron temperature drop inside the sheath tends to diminish the
kinetic correction increasing $\hat{\lambda}$. 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 \in [15, 30] \text{ 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.

Acknowledgement

This work has been supported by the French National Research Agency grant SISTEM (ANR-
19-CE46-0005-03) and has received funding from the Excellence Initiative of Aix-Marseille
University - A*Midex, a French “Investissements d’Avenir” program AMX-19-IET-013. This
work has been carried out within the framework of the EUROfusion Consortium, funded by
the European Union via the Euratom Research and Training Programme (Grant Agreement
No 101052200 - EUROfusion). Views and opinions expressed are however those of the authors
only and do not necessarily reflect those of the European Union or the European Commission.
Neither the European Union nor the European Commission can be held responsible for them.

References

[1] V Philipps, J Roth, and A Loarte. Key issues in plasma-wall interactions for iter: a european approach. *Plasma physics and controlled fusion*, 45(12A):A17, 2003.

[2] Th Pütterich, R Neu, R Dux, AD Whiteford, MG O'Mullane, HP Summers, et al. Calculation and experimental test of the cooling factor of tungsten. *Nuclear Fusion*, 50(2):025012, 2010.
[3] V Ostuni, J Morales, J-F Artaud, C Bourdelle, P Manas, N Fedorczak, R Dumont, M Goniche, P Maget, Y Peysson, et al. Core radiative collapse characterisation and integrated modelling in west plasmas. *Nuclear Fusion*, 62(10):106034, 2022.

[4] S Brezinsek, A Pospieszczyk, G Sergienko, R Dux, M Cavedon, M Faitsch, and K Krieger. Chemically assisted physical sputtering of tungsten: Identification via the 6π→ 6σ+ transition of wd in textor and asdex upgrade plasmas. *Nuclear Materials and Energy*, 18:50-55, 2019.

[5] GJ Van Rooij, Olivier Meyer, S Brezinsek, C Desgranges, D Douai, S Ertmer, A Gallo, C Gil, J Gunn, T Loarer, et al. Tungsten divertor sources in west related to impurity inventory and local plasma conditions. *Physica Scripta*, 2020(T171):014060, 2020.

[6] G Fussman, W Engelhardt, D Naujoks, K Asmussen, S Deschka, A Field, et al. *Proc. 15th Int. Symp. Plasma Physics and Controlled Nuclear Fusion Research*, volume 2, pages 143-1995. IAEA, 1995.

[7] AV Chankin, DP Coster, and R Dux. Monte carlo simulations of tungsten redeposition at the divertor target. *Plasma Physics and Controlled Fusion*, 56(2):025003, 2014.

[8] D Tskhakaya, Mathias Groth, et al. Modelling of tungsten re-deposition coefficient. *Journal of Nuclear Materials*, 463:624-628, 2015.

[9] J Guterl, I Bykov, R Ding, and P Snyder. On the prediction and monitoring of tungsten prompt redeposition in tokamak divertors. *Nuclear Materials and Energy*, 27:100948, 2021.

[10] A Kirschner, V Philipps, J Winter, and U Kögler. Simulation of the plasma-wall interaction in a tokamak with the monte carlo code ero-textor. *Nuclear Fusion*, 40(5):989, 2000.

[11] J Romazanov, S Brezinsek, D Borodin, M Groth, S Wiesen, A Kirschner, A Huber, A Widowson, Markus Airila, A Eksaeva, et al. Beryllium global erosion and deposition at jet-ilw simulated with ero2. 0. *Nuclear Materials and Energy*, 18:331-338, 2019.

[12] HA Kumpulainen, M Groth, G Corrigan, D Harting, F Koechl, AE Jaervinen, B LoMANowski, AG Meigs, M Sertoli, and JET Contributors. Validation of edge2d-eirene and divimp for w sol transport in jet. *Nuclear Materials and Energy*, 25:100866, 2020.

[13] K Shimizu, H Kubo, T Takizuka, M Azumi, M Shimada, S Tsuji, N Hosogane, T Sugie, A Sakasai, N Asakura, et al. Impurity transport modelling and simulation analysis of impurity behavior in jt-60u. *Journal of nuclear materials*, 220:410-414, 1995.

[14] Timothy Reed Younkin, David L Green, Adam B Simpson, and BD Wirth. Gitr: An accelerated global scale particle tracking code for wall material erosion and redistribution in fusion relevant plasma-material interactions. *Computer Physics Communications*, 264:107885, 2021.
[15] Michael Warwick Thompson. Ii. the energy spectrum of ejected atoms during the high energy sputtering of gold. *Philosophical Magazine*, 18(152):377-414, 1968.

[16] A Eksaeva, E Marenkov, D Borodin, A Kreter, M Reinhart, A Kirschner, J Romazanov, A Terra, S Brezinsek, and K Nordlund. Ero modelling of tungsten erosion in the linear plasma device psi-2. *Nuclear materials and energy*, 12:253-260, 2017.

[17] G Falcone. Sputtering theory. *Rivista del Nuovo Cimento della Societa Italiana di Fisica*, 13(1):1-52, 1990.

[18] A Kirschner, D Matveev, D Borodin, Markus Airila, S Brezinsek, Mathias Groth, S Wiesen, A Widdowson, J Beal, HG Esser, et al. Modelling of the material transport and layer formation in the divertor of jet: comparison of iter-like wall with full carbon wall conditions. *Journal of nuclear materials*, 463:116-122, 2015.

[19] D Borodin, S Brezinsek, J Miettunen, M Stamp, A Kirschner, C Björkas, M Groth, S Marsen, C Silva, SW Lisgo, et al. Determination of be sputtering yields from spectroscopic observations at the jet iter-like wall based on three-dimensional ero modelling. *Physica Scripta*, 2014(T159):014057, 2014.

[20] Nicolas Mellet, JP Gunn, B Pégourié, Y Marandet, Céline Martin, and Pascale Roubin. Tungsten erosion by impurities and redeposition: focus on the magnetised sheath. *Plasma Physics and Controlled Fusion*, 59(3):035006, 2017.

[21] D Tskhakaya. One-dimensional plasma sheath model in front of the divertor plates. *Plasma Physics and Controlled Fusion*, 59(11):114001, 2017.

[22] A Kirschner, D Tskhakaya, G Kawamura, D Borodin, S Brezinsek, R Ding, Ch Linsmeier, and J Romazanov. Modelling of impurity transport and plasma-wall interaction in fusion devices with the ero code: basics of the code and examples of application. *Contributions to Plasma Physics*, 56(6-8):622-627, 2016.

[23] A Kirschner, D Tskhakaya, S Brezinsek, D Borodin, J Romazanov, R Ding, A Eksaeva, and Ch Linsmeier. Modelling of plasma-wall interaction and impurity transport in fusion devices and prompt deposition of tungsten as application. *Plasma physics and controlled fusion*, 60(1):014041, 2017.

[24] JP Gunn, S Carpentier-Chouchana, F Escourbiac, T Hirai, S Panayotis, RA Pitts, Y Corre, R Dejarnac, M Firdaouss, M Kočan, et al. Surface heat loads on the iter divertor vertical targets. *Nuclear Fusion*, 57(4):046025, 2017.

[25] I Borodkina, D Borodin, A Kirschner, IV Tsvetkov, VA Kurnaev, M Komm, R Dejarnac, and JET Contributors. An analytical expression for the electric field and particle tracing in modelling of be erosion experiments at the jet iter-like wall. *Contributions to Plasma Physics*, 56(6-8):640-645, 2016.
[26] Jeffrey N Brooks. Near-surface sputtered particle transport for an oblique incidence magnetic field plasma. *Physics of Fluids B: Plasma Physics*, 2(8):1858–1863, 1990.

[27] PC Stangeby. The chodura sheath for angles of a few degrees between the magnetic field and the surface of divertor targets and limiters. *Nuclear Fusion*, 52(8):083012, 2012.

[28] HP Summers. The adas user manual, 2004.

[29] JP Gunn. Kinetic investigation of a collisionless scrape-off layer with a source of poloidal momentum. *Journal of nuclear materials*, 337:310–314, 2005.

[30] Peter C Stangeby. *The plasma boundary of magnetic fusion devices*. Institute of Physics Pub. Philadelphia, Pennsylvania, 2000.

[31] JP Gunn, V Fuchs, and M Kočan. A kinetic model of retarding field analyser measurements in strongly magnetized, flowing, collisional plasmas. *Plasma Physics and Controlled Fusion*, 55(4):045012, 2013.

[32] JP Boris. Relativistic plasma simulation-optimization of a hybrid code, paper presented at fourth conference on numerical simulation of plasmas, naval res. *Lab., Washington, DC*, 1970.

[33] Michael A Lieberman and Alan J Lichtenberg. *Principles of plasma discharges and materials processing*. John Wiley & Sons, 2005.

[34] J Guterl, WR Wampler, D Rudakov, T Abrams, HQ Wang, AG McLean, and P Snyder. Reduced model of high-z impurity redeposition and erosion in tokamak divertor and its application to diii-d experiments. *Plasma Physics and Controlled Fusion*, 61(12):125015, 2019.

[35] C Garcia-Rosales, W Eckstein, and J Roth. Revised formulae for sputtering data. *Journal of nuclear materials*, 218(1):8–17, 1995.

[36] VA Abramov, Yu L Igitkhanov, VI Pistunovich, and VA Pozharov. First wall and divertor plate sputtering in a tokamak reactor. *Journal of Nuclear Materials*, 162:462–466, 1989.

[37] Yonggang Li, Yang Yang, Michael P Short, Zejun Ding, Zhi Zeng, and Ju Li. Ion radiation albedo effect: influence of surface roughness on ion implantation and sputtering of materials. *Nuclear Fusion*, 57(1):016038, 2016.

[38] Rainer Behrisch and Wolfgang Eckstein. *Sputtering by particle bombardment: experiments and computer calculations from threshold to MeV energies*, volume 110. Springer Science & Business Media, 2007.

[39] CA Johnson, SD Loch, and DA Ennis. The effect of metastable atomic levels in neutral tungsten for erosion and impurity transport applications. *Plasma Physics and Controlled Fusion*, 62(12):125017, 2020.