The kinetic theory of prompt redeposition in the case of thin Debye sheath

Understanding the physics of material erosion, migration, and re-deposition on plasma-facing materials is important for the design of fusion tokamak-reactors.^1,2 Prompt redeposition is a process in which a neutral atom sputtered from a divertor target or the first wall of a tokamak returns on the surface right after ionization, before experiencing any collisions with the background plasma. The magnetic field lines in tokamak configuration have a grazing incidence angle on divertor targets. If a sputtered neutral is ionized close enough to the surface (approximately twice the Larmor radius), it can immediately return to the surface due to Larmor gyration in almost parallel to the surface magnetic field. Efficiency of prompt redeposition can be characterized by a prompt redeposition coefficient R defined as the ratio of the redeposited flux to the flux of eroded particles.

Prompt redeposition can be relevant for tokamak plasma-facing materials. In particular, there is an experimental evidence that prompt redeposition is responsible for a large (at least 30%–40%) difference between gross and net erosion of tungsten (W) divertor targets.^3–5 

The Fussmann's equation was tested many times by comparing it to numerical Monte Carlo (MC) simulations^6–8 and seems to be a reasonable first approximation for calculating R. The Fussmann's approach does not take into account influence of the magnetic pre-sheath and Debye sheath on the ions movement, possible multiple ionizations, and other effects, which can be important in divertor conditions.^9,10 Despite that, it is still often used to estimate the prompt redeposition efficiency, e.g., in SOLPS modeling^11,12 because coupling of 2D transport codes to MC simulations of the prompt redeposition can be difficult (although such models do exist, e.g., see Ref. 13). Indeed, from the point of view of 2D plasma codes, only particles, which were not promptly redeposited, matter as they end up in the plasma. Simulations show that R values can be quite large, R > 0.9, meaning that only a small fracture of test particles, 10%, will be relevant in the MC simulations, while the remaining 90% will return to the surface and waste the simulation time. Therefore, further development of approximate analytical expressions for the prompt redeposition is important.

In this work, a new and more universal method is proposed to estimate the prompt redeposition coefficient. We solve the kinetic equation analytically to obtain an analytical expression for the particle distribution function, which then is used to define not only R but also other parameters, such as the average angle of incidence, angular and velocity distributions of redeposited particles, etc. Our expression for R explicitly depends on the angular and velocity distributions of eroded particles and, therefore, is more universal than the Fussmann's result, which is valid for the cosine distribution. Although our consideration at the present moment does not include the sheath electric field and other effects, which are proven to be important by MC codes, we suppose that it is a useful improvement of Eqs. (1) and (2). We also notice that the sheath effects can be included within the framework of our model, though it would make the final equations much more cumbersome.

The MC simulations of prompt redeposition^6,9,10,13 dealt with the cosine distribution of sputtered particles, as anticipated for physical sputtering of pure W divertor targets. However, the angular and velocity distributions of particles sputtered by low energy projectiles near the sputtering threshold are not cosine.^14–16 More complicated distributions are also observed in sputtering of surface structures like fuzz.¹⁷ Using the developed model, we study, for the first time, the influence of the angular and energy distributions of sputtered particles on the prompt redeposition efficiency and parameters of redeposited particles.

Let us consider equations for the distribution functions of neutrals sputtered from a target surface and ions originated from ionization of those neutrals. We suppose that the magnetic field $B = B e z$ is parallel to the target surface, which coincides with the xz plane. The y axis is directed normal to the surface. Both the distribution functions of neutrals f₀ and ions f depend only on y.

In the tokamak divertor configuration, the magnetic field is slightly inclined to the surface. Taking this into account would add f dependence on the coordinates along the surface and make us consider the surface boundaries, which is clearly impossible in our relatively simple analytical approach. If we neglect this additional dependence, one can simply use the parallel (along the surface) component of B instead of total B. Indeed, the Lorenz force component associated with the perpendicular to the surface B component is directed along the surface and does not enter the kinetic equation (4), if there is no f dependence on the coordinates along the surface.

Equation (9) describes a circle of Larmor radius $ρ$ in the xy plane traversed clockwise. The characteristic curve should start at a point $( φ 0 , y 0 )$ where the value of the distribution function, $f ¯$⁠, is known and end at a point $( φ , y )$⁠, where we want to know f. If we are interested in the ions backflux to the surface, we need to consider only the characteristics crossing the surface y = 0 as shown in Fig. 1. Then, a natural choice of $φ 0$ is $φ 0 = 2 π − φ$ as there are no ions leaving the surface; therefore, we simply have $f ¯ = 0$⁠. Integrating (9) and substituting in (11) yields

This result can be intuitively understood as follows. The value of $f ( φ )$ is proportional to the total number of ions crossing the surface at angles $π < φ < 2 π$⁠, see Fig. 1. For a chosen value of $v ⊥$⁠, all such ions come to the surface only along an arc of a circle of radius $ρ$⁠. The circle center position is univocally defined by $ρ$ and $φ$⁠. If $φ > 3 π / 2$⁠, the center of the circle is higher than the surface (left circle in Fig. 1) and is lower for the opposite case, $φ < 3 π / 2$ (right circle in Fig. 1), as the ions gyrate clockwise in our geometry. To find out how many ions move along the circle, we have to integrate the ion source $f 0 / τ$ along it. The number of ions originating in 1 s at a given position y with a given velocity direction $φ$ is $f 0 ( y , φ ) / τ$⁠, and the time to rotate by an angle of $d θ$ is $d θ / Ω$⁠. The total number of ions is obtained integrating along the corresponding arc, and we come to Eq. (13).

The redeposition coefficient is then obtained as $R = Γ / Γ n$⁠, where $Γ$ and $Γ n$ are given by Eqs. (5) and (6).

It is typically assumed that angular and velocity distributions are independent, $F 0 ( v ) = F 0 v ( v ) F 0 φ ( φ )$⁠. In this case, the velocity part of the distribution cancels out, and the redeposition coefficient depends only on the angular part. However, this is the case only if we neglect the influence of the sheath electric field on the ion motion, which could change the distribution function and introduce the dependence on $F 0 v$ in the redeposition coefficient.

If $Ω τ ≫ 1$⁠, the neutrals are ionized much further from the surface than the ion Larmor radius, and $R ≈ 1 / ( Ω τ cos ψ 0 )$ tends to 0. In the opposite case, $Ω τ ≪ 1$⁠, the ionization rate is large, and most of the neutrals are ionized quickly, hence giving R = 1.

This result can be obtained simply as follows. The neutral flux decreases exponentially with y, $Γ = Γ 0 exp ( − y / ( v τ ) )$⁠. All the neutrals which are ionized at $y < ρ$ return to the surface. Their number is $Γ 0 − Γ = Γ 0 − Γ 0 exp ( − ρ / ( v τ ) )$⁠, and we come to Eq. (18).

If $ψ 0 = 0$⁠, Eq. (17) gives $R = 1 / 2$⁠. Indeed, half of the sputtered neutrals fly in the negative x direction (i.e., to the left in Fig. 1), the resulting ions gyrate clockwise, and their trajectories are tangent to the surface, never crossing it, yielding no redeposition flux. The other half goes to the right and crosses the surface immediately after ionization.

These results are illustrated in Fig. 2 where $R ( Ω τ )$ dependencies calculated according to (17) are shown for different $ψ 0$⁠. Notice that if $0 < ψ 0 < π / 4$⁠, the deviation from the perpendicular redeposition coefficient, $R ⊥$⁠, remains less than 10%, and only for larger angles, R becomes sufficiently different, especially at large $Ω τ$⁠.

Calculation of the prompt redeposition for a realistic angular and energy distribution requires numerical integration. In this section, we consider a cosine and “butterfly-like” angular distributions typical for lower energy projectiles. Usually, the angular and energy distributions of sputtered particles are expressed in terms of solid angle and projectile kinetic energy, $F 0 = F E F 0 ψ d E d Ω ψ$⁠, where $d Ω ψ = 2 π sin ψ d ψ$ and angle $ψ$ is referenced relative to the surface normal. This distribution has to be converted in our cylindrical coordinates in velocity space. We have $F E F 0 ψ m 2 d v 2 d Ω ψ = m v F E F 0 ψ v 2 d v d Ω ψ = m v F E F 0 ψ d 3 v$ and $ψ = | π / 2 − φ |$⁠. For example, for the cosine angular distribution, $F 0 ψ = cos ψ$ and $F 0 ( φ , v ) = m v F E ( E ) sin φ$⁠. As was noted before, since integration over angles is separated from integration over velocities, the energy part of the distribution cancels out and does not affect the prompt redeposition coefficient, at least if the sheath influence can be neglected.

Thus, the Fussmann's equations give similar R values for cosine distribution to our result, despite formally the equations look very different. Unfortunately, no detailed derivation of the Fussmann's equations has been published, it is stated only that f(p) is obtained as a ratio of volumes in the velocity space for a given ionization length and cosine distribution of eroded particles, and the integration in Eq. (1) accounts for ionization length distribution.^6,8 However, the ratio of volumes in the velocity space gives the ratio of densities and not the fluxes, which is necessary to define R! Most probably, good agreement of the Fussmann's equations with MC simulations and with our results is due to Eq. (1), where strong exponential dependence on x and p makes the precise form of f(p) to be not that important.

Notice that both curves noticeably differ from the perpendicular distribution coefficient, $R ⊥$⁠, although the difference is only about 10%.

Our consideration ignores the influence of the sheath electric field E on the ion movement. Notice that, in principle, the sheath field can be included in our model. This would only affect the characteristic curve $y ( θ )$⁠. Namely, it would be defined by the solution of a nonlinear differential equation describing movement in crossed E and B fields instead of a simple circle. At the same time, there is no principle difficulty in doing so, and we still could write an expression for the distribution function in a form similar to (13). Interestingly, the effect of the inclined magnetic field is more difficult to take into account within our model because it breaks the symmetry along the surface adding an additional dependence of the distribution function on z (Sec. II).

In the present paper, we do not consider such effects and focus on the case where the most of the sputtered neutrals are ionized outside of the magnetic presheath and neglect the magnetic field inclination. In this case, the sheath does not influence the prompt redeposition coefficient itself, as the electric field does not change the ion flux. However, the electric field still influences the returning ions distribution function. In what follow, we consider its effect under a thin sheath approximation.

The criteria of validity of this approximation can be easily formulated. In the case of an inclined magnetic field, a magnetic presheath occurs, and the electric field E is noticeable at the distances of the order of the main plasma ion gyroradius, $ρ H$⁠, which is larger than the Debye sheath length. The sheath role can be neglected if the sputtered neutrals are ionized outside of the sheath, $λ > ρ H$⁠. This condition limits values of parameter $p ≡ Ω τ$ to $p = λ / ρ > ρ H / ρ$⁠. The energy of the most of physically sputtered atoms equals half of the surface binding energy, E_b. Then, we can estimate $ρ H / ρ = T H M H / ( E b M / 2 )$⁠. Dependencies of $ρ H / ρ$ on T_H for some relevant fusion materials are shown in Fig. 4. The binding energy values are taken from Ref. 18. We see that at low plasma temperatures T_H, the low threshold of $Ω τ$ defined by $ρ H / ρ$ is also reasonably low. These temperatures are typical for the detached divertor but are too low for physical sputtering to occur, except for Li and Be having the sputtering thresholds for H at 5 and 14 eV, correspondingly.¹⁸ At the same time, Li and Sn are suggested to be used in liquid metal divertors where evaporation and thermal sputtering are the main erosion mechanisms with no threshold energies. Redeposited films can also sputter easier than the bulk materials due to weaker bonding between their atoms, compared to better packed crystal lattices of the bulk materials.

It is more convenient to work with the dimensionless variables taking $V 0 = 2 E b / M$ as a natural velocity scale. Then, the sheath strength is characterized by a dimensionless parameter $μ = e U / E b$⁠. As the binding energies for the fusion relevant materials shown in Fig. 4 are of the order of 1–10 eV, and the sheath potential drop ≈3 T_e (for hydrogen plasma), we expect $1 < μ < 100$ for the tokamak divertor conditions.

As an example, let us consider the distribution functions f_s and f₀ for the delta-function angular distribution (14) and Thompson velocity distribution (23), Fig. 5. The neutrals angular distribution is shown in panel (a), the distributions f_s and f₀ are shown in panels (b)–(d) for different $Ω τ$⁠. Parameter $μ = 50$ and $ψ 0 = π / 4$ in all cases.

We see that as $Ω τ$ increases, f₀ peaks blend together. Indeed, at small $Ω τ$⁠, neutrals are instantly ionized, and f₀ is a sum of two delta-functions, with the corresponding angles. Notice that the second peak located at larger $φ$ is slightly higher than the first one as there are more opportunities for an ion to hit the surface at larger $φ$⁠, as one can understand by considering Fig. 1. When ionization weakens, ions can get ionized farer from the surface and may return in a wider range of angles, blurring the peaks of the distribution function.

Therefore, the thin sheath effect on the distribution function is to make ions move closer to the surface normal as one could expect. It causes the peaks to move toward $3 π / 2$ and simultaneously cuts out the ions with too large $φ$⁠, which would correspond to negative expression under square root in (22).

Knowing the distribution function allows, in principle, to calculate any quantity of interest, not only the redeposition coefficient. As an example, we consider in this section the calculation of the average impact angle of ions.

A true value of the prompt redeposition coefficient should be obtained considering the sticking probability of the ions returning to the surface. It is known that the sticking coefficient depends on the angle of incidence. A simple theory¹⁹ predicts that the sticking probability is larger for a grazing incidence, as the particles spend more time near the surface. At the same time, the reverse dependence also exists for some projectiles,²⁰ i.e., the sticking coefficient can either increase or decrease with the incidence angle. Therefore, it is important to calculate the average angle of incidence for the promptly redeposited ions.

As in the previous considerations, we start with the delta-function angular distribution of sputtered particles (14) and the Thompson energy distribution (23), Fig. 6. First of all, one sees that $⟨ φ ⟩ > 3 π / 2$ according to results of Sec. VI, where we have seen that f increases with $φ$ on average. Indeed, for $Ω τ > 1$⁠, there are more trajectories leading to large impact angles. If $Ω τ → 0$⁠, all neutrals are immediately ionized and hit the surface at $⟨ φ ⟩ = 3 π / 2$ due to the symmetry of the initial angular distribution. Evidently, as the sheath potential drop increases, i.e., with $μ$ increases, the ions tend to normal incidence, $⟨ φ ⟩ → 3 π / 2$⁠.

A more realistic case of cosine or low-energy distribution (21) is shown in Fig. 7. Although in both cases, strong sheath influence makes the incidence almost normal, at smaller $μ$⁠, the average angle deviates sufficiently from the normal, by 20° for the cosine distribution. On the contrary, the low energy distribution leads to a very subtle variation of $⟨ φ ⟩$ leaving it equal to normal incidence for all reasonable $μ$⁠. Therefore, in practice, as the variation of the sticking coefficient with $⟨ φ ⟩$ can be strong, we need to account for it, when calculating the prompt redeposition coefficient.

We solve a system of coupled kinetic equations for sputtered neutrals and resulting ions moving in a parallel to the surface magnetic field. By assuming that the sheath is thin, we obtain the distribution function for the ions returning on the surface. Using the obtained distribution functions, we calculate the redeposition coefficient and the average impact angle for various angular distributions of the sputtered particles.

Our results are more general than the Fussmann's formulas for the prompt redeposition (1) and (2). Our equations explicitly depend on the angular and velocity distributions of sputtered particles. It allows the analysis of the impact of these distributions on the redeposition coefficient. We show that R for the cosine and “butterfly”-like distributions typical for low energy sputtering or fuzz sputtering can differ up to 40%, which should be taken into account when considering the prompt redeposition efficiency in a tokamak divertor.

Supposing that the Debye sheath is in thin enough, so that the neutrals are ionized outside of it and the Larmor's radius is larger than the sheath thickness, we have investigated the influence of the thin sheath on the distribution function of the redeposited ions and their average impact angle. We find that the average angle varies significantly with prompt redeposition parameter $Ω τ$⁠, though the variation becomes weaker for larger sheath electric field, as one could expect as the sheath field tends to turn the ions closer to the surface normal.

The main deficiency of our consideration is that we neglect the role of the sheath electric field on the ion movement, except at the end of its trajectory and suppose that the magnetic field is parallel to the surface. These assumptions may not always be valid in divertor conditions. However, we suppose that our results still give useful guidelines for a more detailed investigation of the angular distribution influence on prompt redeposition. We also notice that the MC simulations are not always possible to conduct, for example, for estimation of the prompt redeposition efficiency in 2D transport plasma codes such as SOLPS or UEDGE. In this case, the reported equations can be used, at least as the first approximation.

This work was supported by the Russian Scientific Foundation under Grant No. 24-22-00158 at MEPhI, https://rscf.ru/en/project/24-22-00158/.

The author has no conflicts to disclose.

E. D. Marenkov: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.