Here, we present a numerical fluid plasma model able to capture the enhanced sputtering yield from the Faraday Screen and the Plasma-Facing Components of an Ion Cyclotron Resonance Heating antenna in a fusion machine. The model is a one-dimensional phase-resolved representation of a rectified radio frequency sheath in a magnetic field at an angle with respect to the material surface; the momentum transport of both ions and impurities is computed in the model. The sputtering behavior of the impurities coming off from the wall is obtained from the plasma-material interaction code Fractal-Tridyn. This study analyzes a range of magnetic angles and wave frequencies to parametrically investigate their effect on the energy-angle distributions of the impacting ions and sputtered impurities. Finally, an estimate of the impurity fluxes and of the gross-erosion rate is provided and compared with experimental data available in the literature.

## I. INTRODUCTION

Ion Cyclotron Resonance Heating (ICRH) is envisioned as one of the principal strategies of auxiliary heating power for ITER and future commercial fusion reactors.^{1} The ICRH technique uses radio frequency waves in the ion cyclotron range (typically a few tens of mega-Hertz for magnetic fields of a few Tesla) to heat the plasma by injecting a RF wave into the plasma via an antenna covered by a Faraday Screen (FS) of the conductive material. Proper injection of ICRH power requires the ICRH antenna structure to be placed in close proximity to the plasma edge to attain good coupling between the radio waves and plasma. The close proximity of the biased ICRH antenna combined with the presence of RF waves leads to the emergence of magnetized radio frequency plasma sheaths. These plasma sheaths are mainly driven by radio frequency voltages and are hence known as RF sheaths.^{2–4} The formation of a RF sheath in front of an antenna surface has great implications for the survivability of the Plasma Facing Components (PFC) of the antenna. A radio frequency sheath exhibits rectified voltages much larger than the local plasma temperature, leading to enhanced sputtering of impurities, self-sputtering, hot spot formation, and an overall reduction of the RF coupling efficiency.^{5,6} Furthermore, production of W from the ICRH antenna has been linked to the acceleration of light impurity ions in rectified RF sheaths.^{7,8} Understanding and minimizing impurity sputtering from the antenna surfaces is thus of crucial interest to prevent erosion of ICRH antenna surfaces and edge/core contamination.

Previous experiments on C-MOD^{9,10} and on JET^{11–13} indicate a correlation between localized enhanced impurity sputtering and the presence of RF-sheath potentials on or around active antennas. Plasma surface interaction experiments on the Faraday screen of JET's ICRH antenna coated with beryllium, similar to the Faraday screen which will be used for ITER,^{13–15} showed a significant beryllium influx from optical diagnostics. Erosion experiments on C-MOD using boron wall conditioning on the antenna surfaces estimated a net erosion rate in the range of 15–20 nm/s. Such an estimate would lead to an effective removal of all the boron protective layer in a single 3 MW ICRF discharge operated for only ∼1 s. Despite the differences in size between C-MOD and ITER, the C-MOD ICRH antennas obtain power fluxes of ∼10 MW/m^{2}, in excess of the power fluxes expected on ITER's ICRH antennas. Previous estimates^{16} done assuming a similar erosion rate for beryllium as the boron erosion rate in C-MOD reported that 1000 discharges (400 s discharges) would be sufficient to erode through 1 cm of Be coated on ITER's Faraday screen.

The physical behaviors associated with the emergence of rectified RF sheaths in magnetized conditions^{17,18} and unmagnetized conditions^{19} have been extensively discussed in previous works, with recent experimental investigations now focusing on the subject of impurity sputtering from such systems.^{20–25} Ample evidence and supporting theories^{26} indicate that ICRH antennas enhance the sheath potential on the flux tubes magnetically connected to the antenna surfaces, causing RF sheath potential drops that are significantly greater (of the order of tens to hundreds of T_{e}s) than the classical thermal sheath potential drops (∼3–4 T_{e}, depending on the ion mass). This can be traced back to the misalignment between the current straps of the antenna and the magnetic field. The misalignment leads to the coupling of the slow wave of the ion cyclotron wave with the plasma, generating a nonzero longitudinal electric field,^{5} $E||\u22600$. In the case of ITER, high power ICRH operation (20 MW) demands that the antenna biased voltage should be much larger than the Bohm sheath potential, i.e., $eVRF\u226b3Te$. This bias behaves similarly to a classical thermal sheath (hence, the name rectification), blocking the electrons and accelerating the ions in order to sustain ambipolarity, creating a large rectified direct current (DC) sheath potential in the process. However, in contrast to classical thermal sheaths, ion acceleration is stronger; the ions gain up to hundreds of T_{e}s in kinetic energy, thus potentially causing much more impurity sputtering. Furthermore, the ion dynamics is modulated during each RF cycle. The details of such modulation are extremely relevant for an accurate determination of the ion energy-angle distributions during the RF cycle, and ultimately of the surface response during ion irradiation.

The main goal of the present paper is to provide a one-dimensional numerical model able to account for the sputtered impurities from ICRH plasma-facing components. This has been the subject of microscale^{4,27–31} and whole device^{32–34} modeling attempts in recent years. In the models developed by D'Ippolito *et al.*^{35} and Perkins,^{36} the effects of RF frequency, antenna alignment, and plasma properties were considered, including a preliminary model of the impurity sputtering generated at the ICRH antenna. The model showed an encouraging agreement with experimental data. The same model also predicted the risks of a complete avalanche of self-sputtered screen material impurities^{35} or complete erosion of the target ICRH plasma-facing components.^{36} They also found that coating the ICRH antenna with a low-Z material, such as Be coating, effectively mitigates the probability of self-sputtering avalanche, thanks to the lower self-sputtering yield of low-Z materials. Perkins^{36} also made recommendations on the design variables that would minimize the enhancement of impurity sputtering from the ICRH antenna. However, these models were volumetric (or “zero-dimensional”) approximations of the RF sheath-plasma interaction and failed to capture the relevant physical behaviors occurring inside the RF sheath, necessary to obtain ion energy-angle distributions. More recently, Myra and D'Ippolito,^{5} and Myra,^{4} studied the circuit properties of RF sheaths using an analytical model of higher dimensionality. The Myra and D'Iolito^{5} model is henceforth referred to as the MD model. The MD model developed equations that relate relevant circuit properties to different operational parameters, including wave frequency and plasma properties. However the MD model lacks the ability to estimate the enhanced impurity sputtering yield or the impurity propagation through the sheath.

The main focus of this paper is further extending the MD model, providing a numerical model of the RF sheath that includes impurity generation and transport. The main plasma species in the RF sheath are treated using mass and momentum balances, Boltzmann electrons, and the Poisson equation for the plasma potential. The fluid moments of the ions are then used to sample a population of kinetic ions, treated as computational particles, to obtain ion energy-angle distributions at the wall at the time of impact with the surface. Such distributions are then passed as an input into the sputtering code Fractal-TRIDYN,^{37} for the calculation of the distributions of sputtered impurities. Among the outputs of the code, Fractal-TRIDYN provides the sputtering yields, the reflection yields, and the energy-angle distributions of the particles emitted by the surface during ion-matter interaction. The kinetic distributions provided by Fractal-TRIDYN have been reduced to their zeroth- and first-order moments (impurity density and impurity fluxes) and used as boundary conditions to a set of fluid equations for the impurity species. Thanks to such a model, we have been able to parametrically study a number of features relevant to the ICRH sheaths, namely, the effect of the near-wall plasma parameters (wave frequency, magnetic field angle, and peak-to-peak RF voltage) on the ion impact energy-angle distributions, sputtering yields, and energy-angle distribution of sputtered impurities. Finally, we have obtained estimates of the average flux of sputtered impurities over one RF cycle and the consequent gross-erosion rate.

In Sec. II, we discuss in detail the theoretical model and the assumptions made. The numerical implementation of the model is also discussed. Section III presents the main results of the paper, including the energy-angle distributions of the impacting ions and of the sputtered impurities. In Sec. IV, we discuss the results obtained in Sec. III and calculate relevant erosion rates for comparison with experimental data.

## II. METHODOLOGY

### A. RF Sheath Impurity Model

In this section, we present the equations used to describe the dynamics of the main plasma ions and the impurities released by the wall during a RF cycle. The model is based on a set of one-dimensional phase-resolved fluid equations for the ions and the impurities. Electrons are treated as Boltzmann electrons.

The geometry of the problem is reported in Fig. 1. The domain *x* = [0, *L*] starts in proximity to the wall and is pointed along the normal to the material surface. The origin of the reference frame is placed at the wall (*x* = 0). The domain extends up to the entrance of the magnetic presheath (MPS) (*x* =* L*). The electrostatic field *ϕ*(*x*, *t*) varies only along the *x* axis whereas the constant magnetic field *B* has both *x* and *y* components. The entrance of the MPS (*x* =* L*) is also the point where the ion flow velocity parallel to the magnetic lines is equal to the Bohm acoustic speed *C _{s}*. The simulation domain thus spans from the Debye sheath to the entrance of the magnetic presheath, a region where the ion parallel velocity is always supersonic (larger than the Bohm acoustic speed).

In our model, the 1D equation for local plasma potential with Boltzmann electrons has the form:

where *e* is the fundamental charge, *T _{e}* is the electron temperature in energy units (the Boltzmann constant is implicit), and $ni,e(x,t)$ is the charge density for the ions and electrons, respectively. The electrostatic potential

*ϕ*(

*x*,

*t*) is measured with respect to an unbiased grounded wall (

*x*= 0). In the numerical calculations reported in this paper, we have chosen ions to be deuterium (

*Z*= +1) and impurities to be neutral beryllium atoms, as per ITER ICRH antenna design specifications.

_{i}^{12,15}Equations (1) and (2) are solved with boundary conditions:

where *V _{pp}* is the peak-to-peak amplitude of the oscillating voltage. Sheath rectification is accounted for through an additional condition obtained from the current balances across the sheath, discussed in Sec. II B. The particle densities $nj(x,t)$ and the electrostatic potential

*ϕ*(

*x*,

*t*) are functions of space and time, whereas the electron temperature

*T*is assumed to be constant across the simulation domain (isothermal sheath approximation).

_{e}The continuity and momentum equations for the plasma ions are

where the velocity vector $ui(x,t)=(uxi,uyi,uzi)$ has components expressed in Cartesian coordinates, *ψ* is the magnetic angle (defined as in Fig. 1 as the angle between the magnetic field and the normal to the surface), $\omega ci=eZiB/mi$ is the cyclotron frequency, and *m _{i}* is the ion mass. Equations (1)–(8) form a system of partial differential equations, with boundary conditions given by,

where *n _{L}* is the density at the entrance of the magnetic presheath (

*x*=

*L*) and $|uo|\u2265Cs$ is the magnitude of the ion flow velocity at the same location. Most of our analysis was performed selecting $|uo|$ a little larger than the ion acoustic speed, $uo=\u22121.1Cs$. The tests have revealed that the results are largely insensitive to the actual values chosen, as long as $|uo|$ is strictly larger than

*C*(Bohm-Chodura criterion), so that numerical instabilities are avoided.

_{s}While the simulation is performed in physical units, the following nondimensional units are convenient, especially for comparisons with results from previous RF sheath models:^{5}

where as usual $Cs=Te/mi$ is the ion sound speed, $\lambda D=\u03f50Te/nLe2$ is the Debye length, $\omega pi=nLe2/\u03f50mi$ is the ion plasma frequency, and *E _{th}* is the sputtering threshold energy of the ICRH antenna material. Most of our results will be expressed using such nondimensional units. While $\varphi \u0302$ and $V\u0302pp$ only determine the structure of the Debye sheath, the two nondimensional parameters $\omega \u0302$ (degree of ion mobility) and $\omega \u0302ci$ (degree of magnetization) considerably affect the behavior of the sheath. Their effects will be discussed in Sec. III B.

Once the plasma ions reach the wall, they interact with the material surface. At the energies of interest in an ICRH sheath, the ions gain enough energy to overcome the surface potential barrier and penetrate into the surface lattice. Once inside the lattice, the ions lose energy mainly via Lindhard-Scharff interaction with the lattice electrons, and in part also via large-angle deflections with the lattice nuclei. Such interactions are responsible for a cascade of effects, all accurately accounted for in the Fractal-TRIDYN code.^{37}

In the present work, we have used Fractal-TRIDYN to produce lookup tables of properties relevant to the RF sheath problem. The lookup tables are publicly available as an open-source dataset [Associated dataset available at http://doi.org/10.5281/zenodo.2597591 (Ref. 38)]. The properties of interest are the energy-angle distributions of the sputtered particles, the moments of such distributions, and the sputtering yields $Y(Ei,\theta i)$. In general, the sputtering yield is a function of the energy *E _{i}* and the angle

*θ*of the incoming ions. The total (effective) sputtering yield $Y\xaf$ is found from the weighted integral of the sputtering yield $Y(Ei,\theta i)$ over the normalized impact energy angle distribution function $f\u0302i(Ei,\theta i)$ of the impacting ions

_{i}where the ion distribution $f\u0302i(Ei,\theta i)$ is created assuming a drifting Maxwellian of temperature *T _{i}*, drift velocity

**u**

_{i}and fluid density

*n*. The effective sputtering yield $Y\xaf$ quantifies the gain on a normalized distribution $f\u0302i(Ei,\theta i)$ due to the microscopic yield $Y(Ei,\theta i)$. The flux of sputtered impurities produced by the wall is then given by

_{i}Equation (14) is used as a boundary condition for the sputtered flux at the wall to solve the continuity equation of the impurities

where *n _{I}* is the neutral beryllium impurity density and is initialized to 0 at the start of the simulation. The right hand side of Eq. (15) accounts for sinks of beryllium impurities in the plasma due to electron-impact ionization. All sputtered impurities exit the surface as neutral particles, except for some of the alkali metals of the periodic table, like lithium. In Sec. II B, we describe in detail the boundary conditions for RF sheath rectification. Additional remarks on the sputtering model will be provided in Sec. II C.

### B. RF Sheath Rectification

The total current density arriving at the magnetic presheath entrance can be written as the sum of the ion, electron, and displacement currents

where *J*(*L*, *t*) is the total current density at *x* =* L* and time *t*, *T _{e}* is the electron temperature, $\mu =\nu e/2\pi $ is the effective electron fluid velocity parallel to the magnetic field,

*n*is the normalized density upstream, and

_{L}*ν*is the electron thermal velocity. In order to avoid net DC flow through the system, time symmetry is assumed in our model between the plates i.e.

_{e}Equations (17)–(20) represent the condition that the total current entering the magnetic presheath at phase $\omega t$ must be equal and opposite to the total current at phase $\omega t\u2212\pi $

As a consequence of Eq. (21), the total current integrated over one RF cycle must be equal to zero, expressing the fact that there is no build-up of charge during one cycle. Indicating the time-average over one RF period with angle brackets, we get that the ion flux is

and that the time-averaged displacement current vanishes, $\u27e8\u22022\varphi L\u2202t\u2202x\u27e9=0$. The time-average $\u27e8J(L,t)\u27e9$ of the total current density [Eq. (16)]

where the last term on the left-hand-side is automatically zero under the time average for a periodic wave. Thus, Eq. (23) immediately returns an expression independent of the ion density and the magnetic angle

When *ω* < *ω _{pi}*, the displacement current in Eq. (16) can be neglected without the time average and the total potential (rectified plus harmonics) at the magnetic presheath entrance can be derived

While Eq. (25) is redundant in the presence of Eq. (21), it is numerically enforced during the first half of the RF cycle (*ωt* < *π*) to simulate the initial profile. Equation (21), together with Eq. (25) for the initial *ωt* < *π* phase, is used as a boundary condition on the potential at the magnetic presheath entrance (*x* =* L*) together with the two other conditions presented in Eqs. (3) and (4).

### C. Remarks on the Sputtering Model

An important component of the current model is the inclusion of sputtering yields $Y(Ei,\theta i)$ calculated using the Fractal-Tridyn code^{37} (abbreviated as F-Tridyn). In this section, we provide a few remarks on the procedure used to couple F-Tridyn data in the current RF sheath model.

Each Monte Carlo run of F-Tridyn requires two main types of inputs: (1) atomic and material properties of both the impacting particles and the material surface and (2) the normalized energy-angle distribution $f\u0302i(Ei,\theta i)$ of the plasma ions impacting on the material surface. In order to pass information from the fluid model of the RF sheath to the kinetic model of F-Tridyn, a fluid-to-kinetic conversion is required. We convert the fluid moments (energies and fluxes) produced by the RF sheath model into kinetic distributions using the following approach.

First, we reconstruct a Maxwellian distribution from the two fluid quantities *n _{i}* and

**u**

_{i}, and obtain the ion distribution function in a Cartesian frame $fi(vix,viy,viz)$. Then, we sample a list of particles from this distribution and obtain their kinetic velocities $(vix,viy,viz)$, which can be easily converted into a list of energy and angle coordinates $(Ei,\theta i,\varphi )$, where

*E*is the ion energy,

_{i}*θ*is the ion impact angle with respect to the surface, and

_{i}*ϕ*is the azimuthal angle. The microscopic sputtering yield $Y(Ei,\theta i)$ does not depend on the azimuthal angle, so that coordinate is ignorable

The list of particle coordinates $(Ei,\theta i,\varphi )$ is then binned and normalized, to obtain $f\u0302i(Ei,\theta i)$. Since our major goal is investigating the effect of the energy and angle dependence of the ions on the effective yield, the effective yield was obtained by a weighted integral of the normalized distribution $f\u0302i(Ei,\theta i)$ over the microscopic yield $Y(Ei,\theta i)$

While direct code coupling with F-Tridyn is possible and has been explored in a previous work,^{39} F-Tridyn simulations are computationally expensive at energies larger than > 300 eV. The computational cost and the process of recalculating the sputtering yield for the same set of input parameters make it more efficient to produce a dataset in the form of a lookup table. We produced a dataset covering ion energies from 0 to 1000 eV with 10 eV intervals. We used linear interpolation for data points between intervals to refine the dataset. The associated dataset has been made available at http://doi.org/10.5281/zenodo.2597591 (Ref. 38).

### D. Numerical implementation

The model presented in Secs. II A–II C has been numerically discretized and implemented in a Matlab code. The routines are maintained on a Git repository at the University of Illinois and are freely available upon request. In this section, we briefly describe the numerical methods adopted for the discretization.

The equations of the electric potential, Eqs. (1) and (2), together with the boundary conditions of Eqs. (3), (4), and (25), are solved by finite-differentiation of the Laplace operator and by means of a Newton-Raphson scheme for the nonlinear term deriving from the Boltzmann electrons. The system of equations Eqs. (5)–(8) expressing the ion continuity and momentum is discretized using an explicit upwind scheme and integrated in time until the solution relaxes to a periodic state over several RF cycles. Simulating a small number of RF cycles does not achieve adequate relaxation and significantly overestimates the ion densities at the FS Wall. On the other hand, a large number of RF cycles amplifies the numerical errors leading to the failure of the simulation, while being computationally expensive. The number of grid points necessary for a successful simulation varies depending on the problem parameters *ω*, *ψ*, and *ω _{ci}*. In order to achieve convergence and relaxation to a periodic solution, tuning of the discretization parameters is required, namely: number

*N*of spatial grid points, number

*M*of points per RF cycle, total number of RF cycles, and simulation domain size

*L*.

The procedure used for the implementation of the numerical method was the following. First, we solved the RF sheath model of Sec. II A model using the initial conditions

and initial moments

A typical simulation had to run for at least 10 RF cycles in order to converge to a periodic solution in time.

Good agreement of our implementation with previously published results has been obtained. Figure 2 shows an example of comparison of our model with previous literature results reported in the MD model. Furthermore, the size of the Debye sheath of ∼10*λ _{D}*, defined by the violation of quasineutrality ($ni\u2260ne$), was found to be equal to the values reported in the MD model. Finally, the magnetic presheath was observed to disappear in regimes of weak ion mobility ($\omega \u0302=9$) and magnetic field perpendicular to the wall $(\psi =0\xb0)$, in agreement with previous work. In such conditions, the Bohm-Chodura

^{40}criterion and the classical Bohm sheath criterion are equivalent, and hence, a MPS does not emerge. All the simulations presented in Sec. III are performed in conditions of low ion magnetization parameter, $\omega \u0302ci=0.1$. The effects of ion magnetization $\omega \u0302ci$ have been previously discussed;

^{5}here, we will focus on characterizing the effect of the RF sheath on impurity production.

## III. RESULTS

### A. Overview

The major goal of the analysis presented in this work is characterizing the sputtering behavior of a material wall interfaced with a radio frequency sheath. The computational results presented in this section summarize quantitative simulations aimed at providing a general picture of the surface response under a range of conditions relevant for ICRH and LH radio frequency actuators. Since the sputtering behavior is tightly connected to the macroscopic structure of the plasma sheath, first we outline the behavior of plasma density (both *n _{e}* and

*n*), ion drift velocity, and electric potential as a function of the RF cycle, as obtained from the numerical solution of the model of Sec. II. Then, the calculated energy-angle distributions of the ions impacting on the wall are analyzed in detail, showing the dependence of the ion energy and the ion impact angle as a function of the RF phase. Finally, the energy-angle distributions of the sputtered particles are presented. From their moments, the flux of impurities released by the wall during a RF cycle is obtained and used to estimate the amount of impurities released by the surface per RF cycle.

_{i}The simulations reported here cover a range of magnetic field angles from normal incidence to grazing incidence, with $bx=cos\u2009\psi =1.0,0.6,0.2$ (equivalent to magnetic angles $\psi =0\xb0,\u200953.13\xb0,\u200978.45\xb0$, with *ψ* as defined in Fig. 1), and covering different MPS conditions, going from an absent MPS to MPS larger than the DS. While we simulated a range of peak-to-peak normalized voltages $V\u0302pp=200,100,50$, $V\u0302pp=200$ was used in most simulation cases as it is expected to be relevant for actual fusion operation in ITER.^{12,13} Using an electron temperature *T _{e}* of the order of 3 eV, representative of the typical low-temperature plasma conditions in the scrape-off layer (SOL), the applied RF wall bias voltage is of the order of

*V*= 600 V. While we simulated a wide range of ion mobilities $\omega \u0302=0.07,\u2026,10$, most of our analysis was focused on $\omega \u0302=0.3,1.0,9.0$. These three cases cover significantly different physical regimes and sheath structures, going from highly mobile ions to highly immobile ions.

_{pp}The organization of the sections is the following. In Sec. III B, we report a brief description of the main features of the RF sheath observable from our model. The impact angle distribution and the impact energy of the ions on the FS wall at different times during a RF cycle are presented in Secs. III C and III D, respectively. Section III E discusses the resulting Energy-Angle distributions of the sputtered particles under different physical conditions. Sections III F and III G deal with the first order moment (particle fluxes) of the distribution of sputtered particles.

### B. Sheath Structure

While the structure of a RF sheath is similar to a classical sheath in its constituents, several differences distinguish a RF from a classical sheath.^{5} Figure 3 emphasizes the major features of the sheath structure, which are evident with large $V\u0302pp$. The figure shows a typical solution from the model described in Sec. II. The presence of a rectified potential is the main physical phenomenon that leads to changes in sheath structures. Comparing the average potential upstream $\u27e8\varphi (L)\u27e9\u223cVpp/2$ with the potential in classical sheath cases, $\u27e8\varphi (L)\u27e9\u22483$ V gives a clear indication for the presence of DC voltage rectification. During large parts of the RF cycle, the instantaneous upstream potential $\varphi (L,t)$ [derived in Eq. (25)] is significantly higher than the wall potential $\varphi (0,t)=\u2212Vpp\u2009cos\u2009(\omega t)/2$. As a consequence of sheath rectification, large potential drops appear during periods of high negative wall voltage, as seen in Fig. 3(d). As a result of larger potential drops, the electron density *n _{e}* is stripped away during most of the RF cycle creating a sheath that is almost devoid of electrons. There exist however periods during the RF cycle (centered around $\omega t\u223c\pi $) where the potential drop is significantly smaller, yielding potential drops sufficient to create a more standard thermal sheath.

The presence and size of a MPS are highly affected by the magnetic field angle *ψ* as the MPS has to support a potential gradient that is able to accelerate the ions from $|ux|=cs\u2009cos\u2009\psi $ at the entrance of the MPS to the entrance of the DS given by the Bohm condition $|ux|\u2265cs$. In the cases of a perpendicular magnetic field (*ψ* = 0), the entrance of the MPS and the entrance of the DS coincide ($|ux|\u2265cs\u2009cos\u2009\psi =cs$) leading to the absence of a MPS. The presence of a MPS is detected when a large drop in ion density is visible in the MPS region. The ions are accelerated throughout the DS reaching the FS wall with hypersonic ion normal speeds and higher impact energies.

### C. Ion Impact Energy as a Function of RF Cycle

The ion energy at the time of impact with the wall is a quantity of utmost importance in determining the sputtering response of the surface. Unlike a classical thermal sheath, in a RF sheath the ion energy changes as a function of the RF cycle. Such change in energy has dramatic consequences on the sputtering behavior, since in most of the cases of relevance, the values of the energy fall within a range where the sputtering yield has an exponential variation. This leads to a strongly nonlinear response of the impurity release as a function of the RF cycle.

Figure 4 shows the ion impact energy distribution vs time during one RF cycle at various frequencies $\omega \u0302=0.3,1,3$ and RF voltages $V\u0302pp=200,100,50$. The ion impact energy exhibits strong nonlinear modulations as a function of the RF phase. During the portion of the RF cycle when the wall is negatively biased, the ion peak energy is of the order of the total peak-to-peak voltage plus the thermal kinetic component $eVpp+Ek$. During the remainder of the RF cycle (wall positively biased), the ion energy rapidly drops.

Figure 5 (top) shows the most probable ion impact energy vs time during one RF cycle for various $\omega \u0302=0.07,\u2026,9$ and $Vpp\u0302=200$. Strong nonlinear oscillation in the most probable ion impact energy arise during the RF cycle. The nonlinear oscillations are quite sensitive to the ratio $\omega \u0302=\omega /\omega pi$ (RF frequency over ion plasma frequency). In fact, as the wave frequency *ω* increases with respect to the ion plasma frequency *ω _{pi}*, the effect of the ion inertia becomes evident, effectively reducing the ion mobility across the sheath. At low frequency, the ions have more time to respond to the changing voltage. Higher frequencies lead to a delayed and more flat response of the energy gain. Two main features are evident from Fig. 5 (top). The first is a smoothing of the ion impact energy as a function of time, leading to a decrease in the magnitude of the peak of impact energy. The second is a shift in the ion peak energy toward the center of the RF cycle, toward $\omega t\u223c\pi $. Indeed, ions in higher mobility regimes tend to respond faster to a variation in RF voltage. An increase in $\omega \u0302$ leads to a significant decrease in impact energy during the first half of the RF cycle, with the overall effect of spreading the energy uniformly throughout the RF cycle. In conditions of very low ion mobility ($\omega \u0302>10$), the ion impact energy remains almost constant across the whole RF cycle since the ions are immobile and can only respond to the time-averaged (DC) electric fields.

The nonlinear oscillation of the ion energy during the RF cycle has deep consequences on the sputtering behavior. In Fig. 5 (top), the sputtering threshold of deuterium on beryllium, $Eth=12.6\u2009eV$, is marked with a dashed orange line. Figure 5 (bottom) shows the sputtering yield as a function of time along the same time coordinate, calculated for deuterium on beryllium. Strong nonlinear variations of the yield are clearly evident from the plot. At low frequency, $\omega \u0302=\omega /\omega pi<1$, the ions spend a considerable percentage of the RF cycle near the sputtering threshold. Overall, the increase in the impact energy at the start of the RF cycle does not compensate for the significant decrease in impact energy for the rest of the cycle. The decrease in impact energy below *E _{i}* = 250 eV (close to the maximum of the sputtering curve of deuterium on beryllium) leads to a sharp fall in sputtering yield after the first quarter of the cycle. When averaged over one RF cycle, fewer impurities are sputtered at low frequency. Surprisingly, at high frequency, $\omega \u0302>1$, sputtering is larger than at low frequency. The consequences will be further analyzed in Sec. III F. This can be qualitatively interpreted as follows. While at high frequency the ions do not achieve the same high impact energy of the low frequency case, they spend a larger fraction of the RF cycle above the sputtering threshold, with an overall higher contribution to the yield. The condition $\omega \u0302=\omega /\omega pi\u22481$ (purple curve) represents a regime where the RF frequency approaches the ion plasma frequency. This regime is of most practical interest, as that is the range of frequency that ITER RF launchers will be operating. In such conditions, the behavior is intermediate with respect to the two extreme conditions of high and low frequency. A significant delayed decrease in the sputtering yield occurs during the positive portion of

*V*. Such a decrease is not as shallow as in the low frequency case, but enough to quantitatively decrease the average sputtering yield over one RF cycle.

_{RF}### D. Distribution of Impact Angles

The angular distributions of the impacting ions on the FS wall greatly affect the sputtering yield. Competing forces (electric field force, Hall force, and *E* ×* B* drift) that are a nonlinear function of the RF cycle phase and the operational parameters (magnetic angle *ψ*, normalized wave frequency $\omega \u0302$) dictate the impact ion angular distribution. This creates a strong dependence for the ion impact angle distribution on the RF cycle phase and the operational parameters that could lead to extreme changes in the ion impact angle distribution, from a perpendicular incidence to a back flow of ions in some cases. Such extreme changes in ion impact angle distributions lead to drastic changes in sputtering yield and the energy-angle distribution of sputtered particles. As a consequence, moments of the sputtered particles flux heavily depend on the RF cycle phase and the operational parameters.

Figure 6 shows the normalized impact angle distributions resulting from the relevant physical regimes at different times during a RF cycle. Four cases are shown: (a) magnetic field perpendicular to the surface at low RF frequency ($\psi =0\xb0,\u2009\omega \u0302=0.3$), (b) magnetic field perpendicular to the surface at high frequency ($\psi =0\xb0,\u2009\omega \u0302=3$), (c) magnetic field inclined at an angle with respect to the surface at low frequency ($\psi =78.46\xb0,\u2009\omega \u0302=0.3$), and (d) magnetic field inclined at high frequency ($\psi =78.46\xb0$, $\omega \u0302=3$). Each plot shows four snapshots of the angular distribution at four phases of the radio frequency cycle ($\omega t=0,\pi /2,\pi ,3\pi /2$). The figures show the strong dependence of the ion impact angle distribution on the RF cycle and the operational parameters.

Figure 6(a) (top left) shows the case when the magnetic field is perpendicular to the surface with a low normalized wave frequency ($\psi =0\xb0,\u2009\omega \u0302=0.3$). At the start of the RF cycle $\omega t=0$, the ion impact angle distribution has a sharp bell curve shape around a perpendicular impact angle. The presence of large potential drops creates strong electric fields normal to the wall that dictate the ion trajectory. The electric field forms at a perpendicular angle with the wall; hence, the inclination of the magnetic field determines the role that the *E* ×* B* drift plays. In perpendicular magnetic fields, the *E* ×* B* is identically equal to zero as the magnetic and electric fields are parallel and hence, the magnetic field plays no role. As the potential drop decreases during the RF cycle, the fluid energy of the ions decreases, and becomes comparable to the ion thermal energy. Thermal effects on the impact angle cause a small spread in the distribution between *ωt* = *π* and *ωt* = 3*π*/2. In the high wave frequency case ($\psi =0\xb0,\omega \u0302=3$) [Fig. 6(b,top right)], the ion velocity normal to the wall experiences a smaller decrease in value, remaining well above supersonic speeds. Such a behavior is expected as the oscillations of the ion impact energy during one RF cycle greatly decrease with increasing $\omega \u0302$, a feature that was discussed in Sec. III C.

The steep inclination of the magnetic field gives rise to an *E* ×* B* drift that plays a dramatic role in the trajectory of impacting ions. Figure 6(c, bottom left) show the case of an inclined magnetic field with respect to the surface and a low wave frequency ($\psi =78.46\xb0,\u2009\omega \u0302=0.3$). The presence of an *E* ×* B* drift affects the impact angle distribution in two noticeable ways: shifting the ions away from a perpendicular impact angle and increasing the effect of the Hall force on the ion trajectories. The former change can be explained by the change in th ion path before impact that the *E* ×* B* drift makes. When the magnetic field becomes more oblique, deviation from normal ion impact angles created by the *E* ×* B* drift increases. The latter change can be explained by the increased role that the Hall force plays in more oblique magnetic field cases. The ion guiding center tends to follow a path closer to the magnetic field lines. The change in path creates a deviation in the impact angle distribution center that is directly related to the magnetic field angle. The combination of the physical phenomena listed above creates a heavily amplified effect on the impact angle distribution. At high frequencies and oblique angles, the results shown in Fig. 6(d) are nearly independent of RF phase for reasons already discussed. Compared with Fig. 6(b), the distribution shifts to larger angles corresponding to the change in the magnetic field contact angle.

### E. Energy-Angle Distribution of the Sputtered Particles

The energy-angle distribution of the sputtered particles dictates the kinetics of how impurities leave the surface and flow back into the SOL. Changes in the energy-angle distribution could make a difference between prompt local re-deposition, or contamination of the upstream plasma and potentially of the core. The distributions of sputtered particles obtained from our model can be coupled to external codes handling impurity transport. As a consequence of the dependence of the energy-angle distribution of the sputtered particles on the ion impact angle distribution and ion impact energy, operational parameters (magnetic field angle *ψ*, normalized wave frequency $\omega \u0302$) drastically affect the sputtering response of the surface.

Figure 7 shows the Energy-Angle distribution for different operational parameters at time $\omega t=\pi $. The figure shows the sputtered impurity energy-angle distribution for four different cases: (a) a perpendicular magnetic field at low frequency ($\psi =0\xb0,\u2009\omega \u0302=0.3$), (b) a perpendicular magnetic field at high frequency ($\psi =0\xb0,\u2009\omega \u0302=3$), (c) an inclined magnetic field at low frequency ($\psi =78.46\xb0,\u2009\omega \u0302=0.3$), and (d) an inclined magnetic field at high frequency ($\psi =78.46\xb0,\omega \u0302=3$). In each of the four subfigures, the abscissa is the inclination angle *θ _{s}* of the sputtered particles emitted by the surface, measured as a directional angle with respect to the Y axis. The angle

*θ*is the classical

_{s}*β*angle used in sputtering codes (see, for example, Fig. 1 in Ref. 41 or Fig. 5 in Ref. 37). Note that

*θ*is measured differently than the ion impact angle

_{s}*θ*. The ordinate in each subfigure is the energy at which the particle is sputtered from the wall. The energy-angle distributions of the sputtered particles are a function of the material properties, most noticeably the surface binding energy.

_{i}Sputtered particles leave the surface mostly along oblique trajectories. This is a consequence of the distributions of incoming ions, and the important role that the Hall forces, *E* ×* B* drifts, and sheath E-field play in determining the trajectories of the impacting ions. Figure 7 (bottom left) shows the case with an inclined magnetic field at low frequency ($\psi =78.46\xb0,\u2009\omega \u0302=0.3$). Here, the incoming ions arrive at an angle with respect to the surface, creating a clear shift from normal sputtering and a change in the distribution shape with particles being sputtered preferentially at a sputtering angle of $\theta s\u224860\xb0$. At high frequency ($\psi =78.46\xb0,\u2009\omega \u0302=3$), the ions impact the surface at angles close to $\theta i\u224810\xb0$. As a result, Fig. 7 (bottom right), the distribution of the sputtered particles has a most probable angle slightly deviating with respect to normal incidence.

As the magnetic field angle decreases, the effect of the *E* ×* B* drift significantly decreases. At low frequencies [Fig. 7 (top left)], the particles are sputtered with perpendicular inclinations. When the wave frequency increases and the magnetic field has a perpendicular contact angle with respect to the surface, both the Hall force and the *E* ×* B* drift diminish leading to the electric field force dominating the ion trajectories and causing the ions to impact the surface at near perpendicular angles. In the absence of competing forces, the particles are sputtered preferentially in the perpendicular direction to the surface.

### F. Flux propagation

The amount of impurity fluxes emitted by the RF actuator determines the level of Scrape-Off Layer and Core contamination. In this section, we attempt an estimate of the impurity fluxes in physical units by converting the previous nondimensional model into physical quantities. The calculation is done for a limited number of cases, mainly as a function of frequency and magnetic angles.

Figure 8 shows the effect of the RF modulation of the sputtering yield on the flux of emitted impurities, $\Gamma I(x,t)$. The magnitude of the flux of emitted particles follows the same modulations of the sputtering yield, as determined by the wave frequency $\omega \u0302$. At high wave frequency (case $\omega \u0302=9$ in Fig. 8, top), the sputtered flux does not experience great perturbations during the RF cycle. As a consequence, it remains fairly constant in time. As the wave frequency decreases, the magnitude of the perturbations throughout the RF cycle becomes larger. When the wave frequency *ω* is of the order of the plasma frequency *ω _{pi}* (middle panel of Fig. 8, $\omega \u0302\u22481$), a regime of relevance for ITER operations, the modulation of the sputtering yield causes the magnitude of sputtered flux to experience large oscillations during one RF cycle. Such oscillations are quickly damped out within the first ∼10

*λ*of the sheath from the wall.

_{D}When the wave frequency drops to values lower than the ion plasma frequency (case $\omega \u0302=0.3$ in Fig. 8, bottom), perturbations become even larger in amplitude. In this regime of low frequency, the sputtering yield exhibits highly nonlinear changes during one RF cycle, resembling an on-off behavior. This leads to perturbation of very large amplitude in the flux of emitted impurities.

### G. Average Flux of Sputtered Particles

Over time scales much longer than one RF cycle, the plasma effectively sees a continuous flux of impurities. The time-averaged sputtered flux Γ_{avg} provides a simple scalar quantity useful to compare the net effect of impurity production at the surface of the RF antenna

where $\Gamma I(L,t)$ is the flux of emitted impurities at *x* =* L* as obtained from Eq. (14) and $TRF=2\pi /\omega $ is the RF period. Here, we report an estimate of Γ_{avg}, as obtained from our model, calculated at different frequencies and magnetic angles.

Figure 9 shows the behavior of Γ_{avg} in the low wave frequency range, $\omega \u0302<10$, at three different values of the magnetic field angle *ψ*. As the wave frequency $\omega \u0302$ approaches 0, the flux of sputtered particles Γ_{avg} tends to drop, especially at large magnetic inclinations. Such observation is of easy physical interpretation. At low frequency, the change in RF voltage occurs at a slow pace, such that the ions have high mobility across the sheath potential drop. As the wave frequency increases, the flux Γ_{avg} increases linearly with frequency and decreases with respect to magnetic field inclination. The increase in Γ_{avg} as a function of the frequency is mainly due to an increase in the impurity density *n _{I}* rather than an increase in the impurity drift. Such an increase in the impurity density

*n*with frequency is due to a combination of factors, including changes in the effective sputtering yield, density of the impacting ions

_{I}*n*, and ion impact energy

_{i}*E*. The lowest impurity production is observed at low frequencies and high magnetic inclinations.

_{i}## IV. DISCUSSION

An experimental validation of our model is challenged by the scarcity of experimental data and would require new dedicated experiments. Previous data reported by Bureš *et al.*^{42} were recorded at JET with the new ITER-like wall using beryllium as a plasma-facing component on the Faraday screen. The experiments estimated an average flux of $\Gamma avg=11.9\xd71019m\u22122s\u22121$. Other experiments performed in the US at Alcator C-Mod using an ICRH antenna performed by Wukitch *et al.*^{16} estimated an erosion rate of $15\u201320$ nm/s, using boron as a plasma-facing component on the Faraday screen. Previous impurity sputtering experiments failed to report the operational parameters (e.g., plasma conditions, magnetic field angle *ψ*, etc.) at the Farady screen, which makes any comparison with the modeling predictions if not impossible at least very challenging. A previous modeling effort to validate the impurity sputtering from an ICRH antenna with a beryllium surface was performed and compared to experimental data by Bures^{42} and Wukitch.^{16} Simulations involved a deuterium plasma and a RF wave frequency equal to the plasma frequency (ICRH operational frequency). A range of magnetic field angles was simulated to account for the possible experimental magnetic field angles. Our current model predicts impurity fluxes that are in good agreement with Bureš *et al.*^{42} experimental results, providing values of the average sputtered flux over one RF cycle between $\Gamma avg=9.758\xd71019m\u22122s\u22121$ for a magnetic field perpendicular to the surface ($\psi =0\xb0$), and $\Gamma avg=18.3\xd71019m\u22122s\u22121$ for a magnetic field at large inclination ($\psi \u224880\xb0$). Bureš *et al.*^{42} experiments provided values of $\Gamma avg=11.65\xd71019m\u22122s\u22121$ for beryllium influx from the antennas.

When converted into an erosion rate, $r\u0307=\Gamma avg/NBe$, where *N _{Be}* is the atomic density of beryllium, our calculated value of Γ

_{avg}provides an estimate in the range 80–147 nm/s. The erosion rate estimated by Wukitch

*et al.*

^{16}was calculated for a boron surface, while our simulations all refer to the deuterium/beryllium problem, so the values are not directly comparable.

Projecting our calculated erosion rate to the lifetime of the Faraday screen provides the following considerations. At an erosion rate between 80 and 147 nm/s, a Farady screen covered with a beryllium armor of 1 centimeter of thickness would be consumed in $6.8\xd7104$–$12.5\xd7104$ seconds of continuous operation (about 19–34 h). Such a large erosion rate would completely erode the beryllium armor of the Farady screen in a time not compatible with steady-state reactor operations, putting significant limitations to the operational lifetime of the RF system. Even reducing the operational window of the RF heater would pose significant limitations. For example, for RF discharges lasting 200 s, the antenna surface should be replaced after every $\u223c500$ discharges. While this could be compatible with experimental reactors such as ITER, it would not be acceptable for steady-state operations in a commercial system that needs to operate continuously. Future validation work of the existing model will be necessary in order to address such challenges and identify mitigation strategies.

As a final remark, we would like to point out that the current model does not account for self-sputtering effects on beryllium. This might indeed significantly affect the sputtering behavior in conditions of high re-deposition, due to the lower sputtering threshold and higher yields of beryllium impacting on beryllium. Future works might couple the sputtering model here presented to impurity transport models, in order to obtain a more complete description of the erosion-redeposition processes in front of RF actuators.

## V. CONCLUSIONS

In summary, we discussed a new one-dimensional model based on previous work^{5} describing the physical phenomena of enhanced sputtering from the FS PFC due to the interaction between the PFC and the RF sheaths. The model successfully captured the structure of the sheath, including the DS and the MPS, and their dependence on the different parameters ($\psi ,\omega \u0302,\omega \u0302ci,andV\u0302pp$). The voltage rectification in the sheath, the ion and electron distributions, and the sputtered impurity distributions were also successfully captured. The model has shown that the ion impact energy-angle distribution is heavily influenced by the regime of operation controlled by $\psi ,\omega \u0302,\omega \u0302ci$, and $\u2009V\u0302pp$. In regimes of high magnetic field inclinations $\psi \u224880\xb0$ and low wave frequencies $\omega \u0302=0.3$, the incoming ions achieve highly oblique impact angles at $\omega t=\pi $. The ion distributions shift back toward perpendicular impact angles during periods of high potential drop. As the magnetic field angle *ψ* decreases and the wave frequency $\omega \u0302$ increases, the ion impact angle distribution shifts toward a perpendicular impact angle for larger portions of the RF cycle until reaching a semiconstant impact energy-angle distribution throughout the RF cycle.

Through coupling with F-Tridyn, the effect of the changing RF sheath regimes on the mean flux of sputtered impurities over one RF cycle (Γ_{avg}) is shown in Sec. III G. The mean flux of sputtered impurities Γ_{avg} was found to respond differently to changing wave frequency depending upon the range of wave frequency operation. The erosion rate was calculated for Be FS antennas and was found to be in the range of $\u224880\u2212147$ nm/s. Simulated results have shown good agreement with experimental results, thus providing a basis for possible usage of the code presented here as a reduced-parameter tool in whole-device modeling.

## ACKNOWLEDGMENTS

This material was based upon work supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences, and Office of Advanced Scientific Computing Research through the Scientific Discovery through Advanced Computing (SciDAC) Center for Integrated Simulation of Fusion-Relevant RF Actuators (No. DE-SC0018090-PO97564).