Resonance in radio frequency sheath admittance and enhanced impurity emission near the ion cyclotron frequency

Figure 1 describes the geometry of our model. We simulate the plasma sheath in a 1D domain with a uniform mesh in x, under the justification that the domain is small enough that $\partial_y = \partial_z = 0$ for all equilibrium quantities. While the number of spatial dimensions for the solution of the electric field is kept to one, we maintain three spatial and three velocity components of all kinetic ions. Walls are assumed to be perfectly neutralizing, and so ions which strike the wall leave the simulation. The plasma is maintained by an upstream ion source which maintains a constant volume-averaged density by sourcing ions into the domain with a non-drifting, Maxwellian velocity distribution at some prescribed temperature T_i . New ions are placed randomly in the centre third of the plasma domain. This 'smear' source technique was chosen instead of a spatially uniform ion source to avoid placing new ions directly into the Debye sheath ³ . RF sheaths form at either side of the domain due to time-dependent, periodic, antisymmetric Dirichlet boundary conditions on the electric potential ⁴

$$\begin{align} \Phi(\vec{x} = 0) & = \xi \cos (\omega t) \end{align} \tag{ 1 }$$

$$\begin{align} \Phi(\vec{x} = L) & = \xi \cos (\omega t - \pi). \end{align} \tag{ 2 }$$

The dimensionless quantities used in this work are the following. Both RF and cyclotron frequency are normalized by the ion plasma frequency

$$\begin{align} \hat{\omega} &\equiv \frac{\omega}{\omega_{pi}} \end{align} \tag{ 3 }$$

$$\begin{align} \Omega &\equiv \frac{\omega_{ci}}{\omega_{pi}}. \end{align} \tag{ 4 }$$

This is a natural value to use, as the normalized quantities are typically between 10⁻¹ and 10¹ for ICRH applications and SOL plasmas. $\hat{\omega}$ is referred to as the ion mobility parameter; in the absence of cyclotron resonance, ion sheath admittance maximizes when $\hat{\omega}$ = 1 [26].

Zoom In Zoom Out Reset image size

As usual, T_e is always expressed in electron-volts. Both the peak-to-peak RF voltage $V_\mathrm{pp}$ and electric potential $\phi(\vec{x})$ are normalized by $T_e/e$

$$\begin{align} \Phi(\vec{x}) &\equiv \frac{e\phi(\vec{x})}{T_e} \end{align} \tag{ 5 }$$

$$\begin{align} \xi &\equiv \frac{eV_\mathrm{pp}}{2T_e}. \end{align} \tag{ 6 }$$

This is also a natural normalization for voltage because it is the scaling used for determining the electron density via the Boltzmann relation, meaning both the floating potential of a thermal sheath and the time-dependent rectified potential of a RF sheath scale with $T_e/e$. $\Phi_{\text{bulk}}$ is used to refer to the upstream potential in the middle of the plasma domain. Finally, we normalize various quantities according to the geometry of the magnetic field. If the sheath-producing surface is at x = 0 and $\hat{\mathbf{x}}$ is the direction normal to the surface, then the normal component of magnetic field is defined as

$$\begin{align} \mathbf{b} &\equiv \frac{\mathbf{B}}{||\mathbf{B}||} \end{align} \tag{ 7 }$$

$$\begin{align} b_x &\equiv |\mathbf{b} \cdot \hat{\mathbf{x}}|. \end{align} \tag{ 8 }$$

The angle ψ is then used to describe the angle between the magnetic field and the surface

$$\begin{align} \psi \equiv \arccos b_x . \end{align} \tag{ 9 }$$

2.1.1. Ion sheath admittance definitions.

Ion current density normal to and at the wall for species s is evaluated using the first-order fluid moment of the PIC particle distribution, evaluated on a per-mesh-cell basis,

$$\begin{align} J_{i,s} = Z_s \int q_sf_s(\vec{v}, \vec{x}, t)v_x \mathrm{d}\vec{v} . \end{align} \tag{ 10 }$$

The general expression for the ion contribution to the specific admittance (Ω⁻¹m⁻²) at frequency ω is given by [26],

$$\begin{equation} \text{y}_i = \frac{\langle J_i V\rangle}{\langle V^2 \rangle} - \frac{i \omega \langle J_i \dot{V}\rangle}{\langle \dot{V}^2 \rangle} \equiv G_i + i B_i \end{equation} \tag{ 11 }$$

where V(t) is the voltage of the wall relative to the upstream plasma and $\dot{V} = dV/dt$ and we adopt the $e^{-i\omega t}$ convention for the complex phase. In this section and the rest of the paper, we will refer to the specific ion admittance $\text{y}_i$ as simply the (ion) admittance. The admittance $\text{y}_i$ measures how easily an ion current can flow through the sheath for a given RF voltage. It is conventional to separate real and imaginary components of the RF admittance,

$$\begin{align} G_i & = \mathfrak{Re}(\text{y}_i) = \frac{\langle J_i V\rangle}{\langle V^2 \rangle} \quad &\text{conductance} \end{align} \tag{ 12 }$$

$$\begin{align} B_i & = \mathfrak{Im}(\text{y}_i) = - \frac{\omega \langle J_i \dot{V}\rangle}{\langle \dot{V}^2 \rangle} \quad &\text{susceptance} \end{align} \tag{ 13 }$$

where G_i is the conductance, and B_i is the susceptance of the RF sheath. The polar coordinate $\arctan (B_i/G_i)$ of $\text{y}_i$ in the complex plane is equal to the phase offset between $J_i(\vec{x} = 0, t)$ and $\Phi(\vec{x} = 0, t)$. This phase is largely determined by the ratio of ω to one of the intrinsic ion oscillations such as electrostatic plasma oscillations or cyclotron orbits. When ω is on-par with the dominant intrinsic oscillation, phase-offset is minimized, and $G_i \gg B_i$. Such a sheath is said to be highly conductive for the ions. On the other hand, when $|B_i| \gg G_i$ the frequency of the dominant intrinsic oscillation either exceeds ω (a situation causing ion current to lead RF voltage and $B_i \gt 0$), or is less than ω (a case in which ion current responds to, or lags behind, changes in sheath voltage, leading to $B_i \lt 0$).

In our phase convention, a capacitive susceptance is assigned negative values ($B_i\lt0$), while an inductive susceptance is assigned positive values ($B_i\gt0$),

$$\begin{align} & B_i < 0 \quad \text{Capacitive sheath (ions)} \end{align} \tag{ 14 }$$

$$\begin{align} & B_i > 0 \quad \text{Inductive sheath (ions)} . \end{align} \tag{ 15 }$$

In equation (11) the angle bracket operator $\langle \rangle$ on a generic function of time a(t) denotes the time-average over one RF period $\tau_\mathrm{RF}$,

$$\begin{equation} \langle a(t) \rangle = \frac{1}{\tau_\mathrm{RF}} \int _t ^{t+\tau_\mathrm{RF}} a(t) \mathrm{d}t .\end{equation} \tag{ 16 }$$

For example, the rectified potential in the plasma bulk is obtained from the time-averaged potential over one RF period,

$$\begin{equation} \langle \Phi_{\text{bulk}} (t) \rangle = \frac{1}{\tau_\mathrm{RF}} \int _t ^{t+\tau_\mathrm{RF}} \Phi_{\text{bulk}} (t) \mathrm{d}t .\end{equation} \tag{ 17 }$$

Finally, the admittance $\textit{y}_i$ was normalized with respect to a reference specific ion plasma admittance (Ω⁻¹ m⁻²),

$$\begin{equation} y_{pi} = \frac{\varepsilon_0 \omega_{pi}}{\lambda_D} \end{equation} \tag{ 18 }$$

which yields,

$$\begin{equation} \hat{y}_i = \frac{\textit{y}_i}{y_{pi}} = \frac{\textit{y}_i}{\frac{\varepsilon_0 \omega_{pi}}{\lambda_D}} \equiv \hat{G}_i + i \hat{B}_i \end{equation} \tag{ 19 }$$

where $\hat{G}_i$ and $\hat{B}_i$ are the normalized conductance and susceptance respectively. It is also common expressing the electrical behaviour of a RF sheath using the RF sheath impedance z_i instead of the admittance,

$$\begin{equation} z_i = \frac{1}{\textit{y}_i} \end{equation} \tag{ 20 }$$

which can be simply obtained by inverting the values of RF sheath admittance.

In this work, all plots showing admittance trends were obtained using one hPIC2 simulation per data point. Each simulation used 10⁸ macroparticles, and was configured using a uniform mesh with 200 elements, 100 electron Debye lengths in length. This means that each mesh element contained, on average, 500 000 macroparticles. Such a high number of particles per cell (PPC) allowed for a very low noise evaluation of J_i via equation (10). Two mesh nodes per Debye length were sufficient to sample the Debye sheath.

A time step of 10⁻¹⁰ s was used. The simulation domain was initialized with $\Phi(\vec{x}) = 0$ everywhere. Boltzmann electrons were employed, with T_e held constant at 10 eV. Special care was required to enforce charge conservation in the presence of displacement current for Boltzmann electrons; see appendix A for a detailed explanation of a fast algorithm. Ions were initialized with a uniform spatial density and 3D non-drifting Maxwellian velocity distribution with temperature T_i which, unless explicitly modified, was kept at 10 eV. Each simulation was then evolved in time until $\langle \Phi_{\text{bulk}} \rangle$ and $\langle J_i \rangle$ were constant, averaged over several RF periods. The ion currents were then extracted from the last 10 RF cycles to account for the formation of harmonics in the rectified potential (harmonics in RF sheath potentials were previously seen in the sheath model of Lieberman [27]). Real and imaginary components of the ion sheath admittance were obtained using equation (11).

Each simulation used 1 MPI process, each with 1 AMD Milan CPU with 64 OpenMPI threads, using 256GB of DDR4 RAM, and 1 Nvidia A100 GPU with 6912 CUDA cores and 40 GB of 1935 GB s⁻¹ GPU memory (312 TFLOPS). The combined processing power of the GPU and OpenMP threads allowed for very high fidelity simulations (500 000 PPC) to simulate to ion equilibrium time in a matter of a few hours. Parameter scans were made much more feasible by scheduling tens of such simulations at a time thanks to the UIUC NCSA Delta Computing Center. The combined data from these experiments was obtained using super computing GPU-accelerated hardware which would require hPIC2 over 10 years of run time in single-core CPU mode.

This section introduces observations and results on the formation of an admittance resonance in a magnetized RF sheath at oblique magnetic field angles. Specifically, a resonance is observed in the ion sheath admittance when the RF is in the proximity of the ion cyclotron frequency. For these simulations, a pure deuterium plasma was used. RF was fixed at $f_\mathrm{RF} = 55.9$ MHz, which is in the range of ICRH frequencies used for example in the WEST tokamak [28]. The magnetic field strength was the scaling parameter for each simulation.

Figure 2 shows the behaviour of the real (figure 2(a)) and imaginary (figure 2(b)) components of the RF sheath ion admittance, and of the rectified potential $\langle \Phi_{\text{bulk}} \rangle$ in the plasma bulk (figure 2(c)). The plot is presented as a function of the normalized frequency $\omega / \omega_{ci}$ (calculated by inverting the magnetic field value for each case and scaling by $m_i \omega / e$), over an interval of frequencies spanning from one-half to three ω_ci 's.

Zoom In Zoom Out Reset image size

The plot includes data from four magnetic field inclinations: 30^∘, 45^∘, 70^∘, and 80^∘. Several outstanding features emerge from the plots. First of all, at small magnetic inclinations (see the $\psi = 30^\circ$ case in figures 2(a)–(c), brown curves) the conductance (figure 2(a)) exhibits a local maximum at the ion cyclotron frequency, $\omega/\omega_{ci} = 1$. In the RF sheath ions are resistive, with peak conductance of the order of $\hat{G}_i \sim 10^{-2}$, that is only 1% of the ion plasma conductance. At those small inclinations, the susceptance (figure 2(b)) remains positive, meaning that the RF sheath exhibits inductive behaviour. However, as magnetic inclination increases, the peak in conductance becomes more pronounced, and the susceptance develops a pronounced well, with a transition from inductive ($B_i\gt0$) to capacitive ($B_i\lt0$) behaviour. At larger inclinations, and for frequencies in proximity to the ion cyclotron frequency, $\omega \sim \omega_{ci}$, the RF sheath becomes a better conductor, with peak conductance on the order of $\hat{G}_i \sim 10^{-1}$, that is 10% of the reference plasma ion conductance.

From the trends observed in figure 2 we can define a resonant frequency ω_r :

$$\begin{align} \omega_r \equiv \{\omega/\omega_c \,| \, \mathfrak{Re}(\hat{y}_i(\omega)) = \max(\mathfrak{Re}(\hat{y}_i))\} \end{align} \tag{ 21 }$$

which is the frequency where the conductance has its maximum. At large angles, a capacitive response ($\mathfrak{Im}(\hat{\textit{y}}_i) \lt 0$) is seen for $\omega \lt \omega_r$ where the effect is not washed out by broadening. Furthermore, a second resonance appears near the second cyclotron harmonic, $\omega = 2\omega_{ci}$ for the highly oblique case, $\psi = 80^\circ$. Unexpectedly, at those large angles a resonance peak develops also at the half-resonance, $\omega = \omega_{ci/2}$. This is a rectification phenomenon, arising from a non-linear interaction between the twice-oscillating floating potential and the cyclotron frequency.

One striking phenomenon that occurs during resonance is sheath depletion, which is characterized by a fraction of the RF period during which the ion density n_i drops an order of magnitude, and J_i drops to zero, at the wall. Sheath depletion occurs when the RF sheath excites ion gyro-motion to the point that the average Larmor radius is a significant fraction of the sheath width. As excitation of gyro-motion scales with $\sin\psi$, this effect is observed at oblique angles.

Figure 3 shows the normal component of ion current density versus time for the semi-oblique case of $\psi = 70^\circ$. In figure 3(a) we see a special low-frequency case of $\omega/\omega_r \sim 0.5$. Here, ions complete two gyro orbits per RF period. The phase-shift closely resembles the current–voltage relationship of a capacitor, with the small fundamental component of the current $+\pi/2$ out-of-phase with V. This also makes sense if we look at figure 2, where for the corresponding frequency, $\mathfrak{Re}(\hat{y}_i) \ll 1$. In figure 3(b), we see the current profile for $\omega/\omega_r \sim 0.86$. Here, $\mathfrak{Re}(\hat{y}_i) \gt 0$ and $\mathfrak{Im}(\hat{y}_i) \lt 0$, meaning the phase-shift between J_i and V is less than $+\pi/2$ yet still greater than 0. Higher harmonics have formed in $J_i(t)$. In addition, a peculiar phenomenon has begun to develop: for a brief stretch of time per cycle, the ion current at the wall is positive, flowing back into the plasma. Increasing ω even further to resonance (figure 3(c)) causes the current to change phase toward a mildly inductive regime, with most of the current amplitude coming in the form of conductance. This frequency, which represents the peak of ion admittance in the resonant regime, results in a phase shift of the ion current such that $\mathfrak{Im}(\hat{y}_i) \sim 0$, meaning the sheath is highly conductive to ions. At the frequency depicted in figure 3(c), the gyro-orbit is maximally driven by the RF voltage and the amplitude of $J_i(t)$ is greatest. Due to finite Larmor radius effects, $J_i(t)$ has departed from a harmonic profile. The current profile has formed a pedestal and a portion of it is positive, meaning the ion drift velocity briefly points from the wall back into the plasma. For $\omega \gt \omega_r$, the phase continues to shift, asymptotically approaching $-\pi/2$. For $\omega/\omega_c \sim 2$ (figure 3(d)), the ion contribution to the sheath is almost purely inductive. Above this frequency (not shown), $J_i(t)$ decreases in amplitude and does not change in phase.

Zoom In Zoom Out Reset image size

Figure 4 demonstrates the deformation and inversion of IV waveforms in phase space at different resonances for $\omega = \omega_r$. Ion currents are normalized by their amplitude

$$\begin{align} J_{i,\text{norm}}(t) \equiv \frac{J_i(t)}{\max(J_i) - \min(J_i)} . \end{align} \tag{ 22 }$$

While this normalization is appropriate for comparing the shape of the IV waveforms across simulations with potentially varying amplitudes of $J_i(t)$, the reader should take note that it has the effect of encoding the magnitude of $J_i(t)$ into the y-offset of $J_{i,\text{norm}}$. Currents whose amplitude is small compared to their magnitude will result in $J_{i,\text{norm}} \lt -1$ for more time steps. On the other hand, currents whose amplitude is a significant fraction of their magnitude will result in $-1 \lt J_{i,\text{norm}} \lt 0$. In other words, the closer $\langle J_{i,\text{norm}} \rangle$ is to 0, the greater the ion admittance is.

Zoom In Zoom Out Reset image size

The purpose of this transformation is to emphasize the profile of $J_i(t)$ in phase space. Figure 4 gives a clear demonstration of the anharmonic effects on ion current across the resonance with increasing RF voltage. Since the RF voltage provides the driving force, sheath depletion becomes more dramatic with increasing amplitude ξ (recall that $\xi \equiv V_\mathrm{pp}/2T_e$ is the normalized RF peak-to-peak voltage). Note that, for the plot of figure 4, ion current over ten RF periods is plotted, thus the thickness of the line is an indication of the variance in $J_i(\varphi)$ where ϕ is a fixed RF phase. As we start from a low RF voltage of ξ = 1.25, $J_i(t)$ is entirely negative throughout one period, flowing toward the wall. Doubling the voltage to ξ = 2.5 increases the amplitude of $J_i(t)$ by shifting the low-current tail of the contour upward. This can be attributed to an increase of both Larmor radius and gyrocentre motion in the RF electric field ⁵ . As some ions miss the wall in their orbit, they are soon propelled back into the plasma, contributing a positive amount to the overall ion current. As we double voltage again to ξ = 5, the contour begins to invert. Now, near the high-current side of the contour, we begin to see a steeper curve. Ions have increased their gyromotion to the point where they are more directly striking the wall. However, the increase in Larmor radii means that the fraction of ions which miss the wall are even more propelled back into the plasma, leading to $J_i(t) \gt 0$ for a short range of time steps. We can see an even more dramatic effect as we increase RF voltage to ξ = 25. Examining the high-current side of the RF period first (the left-hand side of the contour) we see a very sharp drop from $J_i(t) \approx 0$ down to maximum current. The ions have been excited to such a degree that a single Larmor orbit is sufficient to traverse the entire sheath. The ions with the largest radii strike the wall and leave the simulation. Soon after, there remains a fraction of the ion population with a slightly smaller Larmor radius, yet still large compared to the sheath width. These ions curl back out of the sheath, traversing the sheath directly. Enough ions leave the sheath during this period that the sheath becomes entirely devoid of ions, as evidenced by the long flat portion of the contour near $J_i = 0$ for ξ = 25.

The implication of sheath depletion is that it can synchronize the direction of ions in their orbits, significantly increasing ion admittance compared to other cases where the gyro-radius is not a significant fraction of sheath width. However, in order for sheath depletion to occur in the first place, the sheath must be quite narrow. The next section demonstrates how a narrow sheath width actually presents a damping effect to the RF driving force on the ions.

It was observed that, in certain conditions, the resonance peak broadens, and ω_r shifts below 1. This phenomenon, referred to as admittance damping, is controlled by plasma density and ion temperature. These observations can be partly explained by the modification of the sheath width which plays an important role in setting the ion transit time across the sheath, as does the ion temperature.

Figure 5 highlights damping as a function of ion temperature for the semi-oblique case of $\psi = 70^\circ$. T_e was kept constant at 10 eV. As ion temperature increases from 5 eV to 1 keV, the real part of admittance broadens slightly and decreases significantly in amplitude by more than 90%. The imaginary admittance profile decreases by more than one order of magnitude and does not appreciably broaden. The source of the damping in this case lies in the inability of the RF sheath to significantly alter the trajectories of cyclotron resonant ions with large thermal velocities. A resonance will occur when $\omega/\omega_{ci}\sim 1$ regardless of how large T_i is since, at that frequency, particles complete one gyro-orbit per RF period; however, as explored later, the time the ions spend in the sheath is proportional to the sheath width divided by the thermal velocity, making the relative energy gain from the cyclotron resonance small at high T_i . Finally, the admittance peak at the second harmonic ($\omega \sim 2\omega_{ci}$) becomes more pronounced with greater ion temperature. This is due to the fact that power absorption at the second harmonic increases with T_i (the ion acceleration due to ion cyclotron resonance is proportional to a Bessel function $J_m(k_\perp v_\perp / \omega_{ci})$, where m is the harmonic number [29]).

Zoom In Zoom Out Reset image size

Figure 6 demonstrates damping as a function of plasma density $\langle n_e\rangle$. This is one of the most important results of this paper. The effect of increasing plasma density is to shift the resonant frequency ω_r downward and to significantly broaden the resonance peak. For low density, and hence low damping, the conductance is nearly zero away from the narrow resonance. For higher density, an increasing fraction of the overall admittance is real, coming in the form of conductance. Note that the normalized frequency $\hat{\omega}$ changes dramatically due to the change in $n_{i,\text{SE}}$. The rectified potential $\langle \Phi_{\text{bulk}}\rangle$ stays relatively constant throughout this parameterization, meaning observed effects are not the result of a change in sheath voltage.

Zoom In Zoom Out Reset image size

The admittance resonance profile deviates from that of a classically-driven, damped harmonic oscillator. This is because the sinusoidal driving force of the E-field which excites particles' Larmor orbits is superimposed on all the other forces at the sheath entrance. In particular, there is a background admittance resulting from these forces which can be appreciated by examining the nearly unmagnetized limiting case of ψ = 30 in figure 2. The forces in this case combine such that the sheath is inductive to ions, causing a nonzero positive susceptance. The resulting shift in $\mathfrak{Im}(\hat{y}_i)$ affects the admittance resonance profile for all resonant cases in magnetized conditions.

One diagnostic we can employ to better understand the mechanism behind damping is the ion sheath transit time (STT). The single-ion STT is defined as the time it takes a particle to travel from the sheath entrance to the wall at x = 0. The sheath entrance location is defined by the Bohm acoustic speed:

$$\begin{align} \Delta(t) &\equiv \displaystyle\left\{\hspace{2mm} x(t) \hspace{2mm} \middle| \hspace{2mm} u_{x,i}(x,t) = -c_s \hspace{2mm} \right\} . \end{align} \tag{ 23 }$$

Distributions in ion STT at various charge densities can be found in the right-hand column of figure 7. Each row corresponds to one hPIC2 simulation run at the resonant frequency ω_r for one of the densities in figure 6, repeated here: 10¹⁵ m⁻³, 10¹⁶ m⁻³, 10¹⁷ m⁻³, and 10¹⁸ m⁻³. In the left-hand column, we see typical particle trajectories across the sheath where the excitation of the gyro-radius produces a cone-like trajectory. The sheath width narrows with increasing charge density in proportion to the Debye length ($\Delta \propto n_{i,\text{SE}}^{-1/2}$). As the sheath narrows, ions spend less time in the sheath. But the RF driving force exciting ion gyro-orbits is only present in the sheath. Therefore, at higher charge density, ions experience a resonant excitation for a shorter period of time.

Zoom In Zoom Out Reset image size

The broadening of the ion admittance resonance due to STT may be captured by the dimensionless factor $\gamma = \tau_c/\text{STT}$, where τ_c is the gyroperiod, which can be used to generalize a fit-function for ion sheath admittance (developed by Myra [26] equation (54)), substituting $\hat{\omega} \longrightarrow \hat{\omega} + i\gamma$:

$$\begin{align} \textit{y}_i \propto \text{f}_{yi} \equiv ib_x^{3/2}\frac{1}{\hat{\omega} + i\gamma}\frac{(\hat{\omega} + i\gamma)^2 - b_x^2\Omega^2}{(\hat{\omega} + i\gamma)^2 - \Omega^2} \end{align} \tag{ 24 }$$

Plugging in the average experimental STT to obtain γ in equation (24) for each of the four densities in figure 6 successfully reproduces the observed admittance broadening. Figure 8 highlights the damping trend for these four densities (as well as a few additional ones) by comparing the full-width half-max (FWHM) to γ. The trend shows excellent agreement with a linear relationship.

Zoom In Zoom Out Reset image size

The broadening of the peak in ion admittance resonance described in this subsection is sufficiently described by inverse sheath transit time and equation (24). As noted earlier, a shorter time spent in the resonant region causes a resonance broadening. This concept is consistent with a number of other resonant phenomena across multiple fields throughout physics. Similar arguments are made in quantum mechanics to predict the broadening of spectral line widths based on the lifetime of the excited state [30–32], as well as in nuclear physics to calculate the half life of highly unstable isotopes [33, 34].

Our observations of ion admittance resonances sheds light on the physical mechanisms responsible for the ion dynamics within the sheath. The consequences on ion energy enhancement, material sputtering and erosion are considered in the next sections.

IEADs of the ions impacting a material surface are a necessary tool for predicting PMI, particularly sputtering. An IEAD is a 2D distribution function of ions crossing the simulation boundary (in this case, representing a material surface) where the degrees of freedom are ion energy and ion angle Θ from the surface normal. In most cases, kinetic plasma models such as particle-in-cell codes are required to accurately analyse PMI because they can provide a detailed IEAD for each surface. The case of an RF sheath in an oblique field near the ion cyclotron frequency presents an extreme case where the rectified potential, which is usually the standard for estimating material sputtering, does not appreciably change across resonance while the IEAD profile changes significantly.

Figure 9 shows D⁺ IEADs below, across, and above resonance for a high voltage, semi-oblique case. Note here, as in the previous section, the scanned parameter was the B-field, with the RF fixed at $\omega = (2\pi)55.9$ MHz. Here we gain a further understanding of the effects of admittance resonance on PMI. For the case of $\omega/\omega_r = 0.52$, the ion cyclotron frequency is nearly double the RF. Ion current peaks twice per RF period, each time accelerated by a slightly different potential due to the change in RF phase. The result is a bifurcation of the distribution peak (figure 9(a)). As we approach resonance ($\omega/\omega_r = 0.86$), the peak in the distribution shifts slightly higher in energy, while the low-energy tail shifts toward lower Θ, representing a more direct impact at the wall. At resonance (ω = ω_r ), the peak of the ion distribution at the wall decreases to ∼20^∘, while the most probable energy more than doubles up to 60 $\mathrm{eV}/T_e$. Elias et al found that, for a range of incident angles, RF sheaths produce IEADs with a long band ranging from E = 0 to approximately $E = \xi$ ([23], figure 4). For both the $\omega/\omega_r = 0.86$ and $\omega/\omega_r$ = 1 cases in figure 9, this band is visible in the background (the light green shading between 30^∘ and 60^∘, in this case up to 30$T_e/e$). One interesting finding by Elias et al was that this energy band disappears at high frequencies (${\unicode{x2A7E}} \hat{\omega} = 2.0$). However, despite the relatively low ion mobility of $\hat{\omega} = 3.2$ for the cases observed in figure 9, this energy band persists near resonance. This exemplifies the ability of cyclotron resonance in the sheath to overcome any lack of ion mobility from plasma oscillations. The main peak, however, is a structure much higher in energy, up to ∼3ξ for $\omega = \omega_r$.

Zoom In Zoom Out Reset image size

Lastly, in the region at and below resonance (figures 9(a)–(c)), a critical angle $\Theta_c$ develops below which no ions strike the wall. This is peculiar because one would expect resonance, which produces the most energetic ions, to accelerate a greater portion of ions to direct impact angles. The width and vertical extent of this depletion band can be explained by the magnetization parameter, defined as $\lambda_D / \rho_s$, where $\rho_s \equiv \text{c}_s/\omega_{ci}$ is the ion sound radius. The magnetization parameters for each of the four cases in figures 9(a)–(d) are, respectively, 3.36, 2.025, 1.749, and 0.851. For figure 9(a), which is the most magnetized plasma, ions are pulled closest to the direction of the B-field (70^∘) and the depletion band is largest, while for the least-magnetized case in figure 9(d), the perpendicular electric field provided by the sheath dominates and the depletion band vanishes.

The effect of admittance damping on the IEAD is highlighted in figure 10. Here we see the IEADs plotted at the resonant frequency $\omega = \omega_r$ for a range of plasma densities. For all densities, the distribution peak at resonance takes a triangular shape, with a long vertical section at direct angles stretching over a range of energies, as well as an elongated mono-energetic tail spanning a range of angles. When damping is at a maximum ($\langle n_e\rangle$ = 10¹⁸ m⁻³), the IEAD is a composite of two peaks. The first peak is the resonant triangle we see for all other densities at resonance, here significantly shifted downward as a result of damping. Superimposed onto this peak is another, narrower peak of ions which impact the wall at greater inclination. This second peak can be explained by the severe truncation of the gyro-orbit due to the small sheath size ⁶ . In addition, at the highest density case of $\langle n_e\rangle = 10^{18}$ m⁻³, the critical inclination $\Theta_c$, observed for all other resonant IEADs, has disappeared, and ions are once again able to directly strike the wall ($\Theta = 0^{\circ}$). This can once again be explained by the ion magnetization parameter, which is at a minimum for the high-density case of figure 10(d), where the magnetization is weakest, enhancing the effect of the sheath electric field accelerating ions directly into the wall over the effect of the ion motion parallel to B. Ultimately, damping has the effect of narrowing the peak of the distribution in energy-angle space.

Zoom In Zoom Out Reset image size

The doubling of ion impact energy at the wall throughout resonance (figure 9) calls for further investigation into material sputtering. Specifically, we wish to investigate the effect of resonance on material sputtering inside a tokamak. The question is where precisely such resonant RF sheaths would form in a tokamak. The formation and coupling mechanisms of RF sheaths on hardware far from ICRH actuators is well known [3]. Assuming a fixed ICRH frequency, resonance would occur at different locations depending on the species mass and charge, as well as the magnetic field strength at one of these sheaths. In a fusion plasma, there is a wide range of species, including oxygen impurities as well as trace amounts of neon from gas puffing near the divertor. Due to the $1/R$ scaling of the toroidal B field in a tokamak, resonances will therefore occur in far-field RF sheaths at various major radii. The only exception to this broad generalization is the set of locations with large plasma densities, since we have shown in section 3.3 that charge density damps admittance resonance. In divertor plasmas, this eliminates locations along the divertor, specifically the near the strike points (the plasma density in ITER is predicted to be as large as $2.5 \times 10^{21}$ m⁻³ [35] at the inner and outer strike points). Enhanced sputtering from cyclotron resonance is therefore not expected in the private flux region or near SOL and we may limit our investigation to far SOL surfaces.

Motivated by these considerations, we ran a series of multi-species RF sheath simulations with hPIC2 [25], and passed the resulting IEADs to RustBCA [24], a binary collision approximation code capable of modelling ion–material interactions such as sputtering, reflection, and implantation. Simulations were run for an isothermal plasma, $T_e = T_i = 10$ eV, at low-density, $n_e = 10^{17}$ m⁻³, with atomic composition approximately representative the far SOL and first wall in ITER-like conditions [36], namely, deuterium–tritium in equimolar concentration, $n_D = n_T = 0.443n_e$, helium-4 from fusion ash (${\lt}1$%) in both charged states ⁴He⁺ and ⁴He²⁺, helium-3 (³He²⁺) at 3% concentration used for minority heating [37], as well as other trace impurities (${\lt}1$%) such as oxygen, neon, and beryllium including all relevant charge states. The complete set of ion densities used in these simulations can be found in table B3. The magnetic field was left as a parameter and scaled from 1 T to 12 T to study the resonance behaviour, with magnetic field inclination $\psi = 70^\circ$ with respect to the surface, sufficiently oblique to induce resonance and demonstrate its effects on sputtering. The RF wave in proximity to the surface was assumed having peak-to-peak amplitude V$_\mathrm{pp}$ = 50 V (ξ = 2.5) and frequency of f = 54 MHz. Two walls were considered, made of pure tungsten and pure beryllium. A ZBL inter-atomic potential [38] was used in RustBCA for ion surface interactions [39]. Note that such a purely repulsive potential ignores surface chemistry effects, which may increase overall sputtering. The complete set of material properties for beryllium and tungsten target materials are reported in tables B1 and B2, respectively.

Figure 11 illustrates the calculated sputtering yields resulting from the hPIC2-RustBCA simulations, for the case of tungsten wall. The different panels of the figure are divided by impacting ion species. Note that here and in the following sputtering yields are denoted by Y_i as is conventional, which should not be confused with the RF sheath ion admittance $\textit{y}_i$. In figure 11, peaks in Y_i can be seen at the resonant magnetic field strength for each species. Two results are particularly striking. The first is the behaviour of low-mass, single-or-doubly charged species such as deuterium (figure 11(f)), ³He²⁺ (figure 11(e)), and ⁴He²⁺ (figure 11(d)). Y_i increases significantly at resonance from these incident ions, more than tripling in the case of deuterium (resonant B-field strengths are indicated by solid circles). While the overall sputtering yield of tungsten from incident deuterium is low (${\unicode{x2A7D}} 10^{-3}$), it must be taken into account that the flux of deuterium at the divertor is the greatest of all species. Such a significant increase of Y_i from deuterium has implications for gross erosion and total impurity production. The second effect is the behaviour of medium-mass, highly-ionized species such as Ne¹⁰⁺ (figure 11(a)) and O⁸⁺ (figure 11(b)). For these ions, a local minimum in Y_i is present just below resonance, most dramatic for high charge states such as Ne¹⁰⁺ and O⁸⁺. This local minimum results from the off-resonance magnetization of species crossing the sheath, causing them to follow the oblique magnetic field lines, thus impacting the tungsten wall at greater inclination where the sputtering yield is lower. Moreover, for these highly ionized species, resonance increases Y_i only slightly. This is due to the fact that the baseline Y_i is already quite large for such highly-energetic ions. Note that the apparent broadening of the resonant peaks of Y_i with decreasing charge state of neon, oxygen, and beryllium is not a result of increased damping; rather it is simply the result of scaling with the magnetic field instead of $\omega/\omega_{ci}$.

Zoom In Zoom Out Reset image size

Figure 12 reports the sputtering yield for the case of a beryllium wall. The resonant response of beryllium Y_i is more or less the inverse of the tungsten sputtering yields, with local minima occurring at resonance rather than peaks. This can be explained by the general rule of thumb that, when the incident ion is more massive than the material atom, the sputtering yield is angle-dominated as opposed to energy-dominated, and peak sputtering yield occurs at oblique angles of incidence [22]. This is why beryllium Y_i from deuterium (figure 12(e)) only has a slight decrease but neon (figure 11(a)) has a significant decrease at resonance.

Zoom In Zoom Out Reset image size

This phenomenon may be more concretely appreciated by examining figure 13, which compares the energy-and-angle dependence of beryllium with tungsten, in this case for incident oxygen. From these results one clearly sees that the sputtering yield for oxygen on beryllium is significantly more dominated by incidence angle than is the case for tungsten. Sweeping mono-energetic 400 eV oxygen, for example, across angles of incidence from 0^∘ to 60^∘, would only slightly vary the sputtering yield of tungsten, while for beryllium the sputtering yield would peak at 70^∘ and decrease up by more than a factor of 5 as the angle of incidence approached the normal. Overlayed on the contour plots in figure 13 are the distributions of O⁸⁺ at the wall from the hPIC2 far SOL RF simulations. One sees that the effect of resonance is to increase the energy slightly while bringing the impact distribution closer to normal incidence. These variations, along with the sputtering yield contours, explain the difference in resonant Y_i response between tungsten and beryllium in figures 11 and 12.

Zoom In Zoom Out Reset image size

The tungsten sputtering results (figure 11) indicate that several important species produce peaks or significant increases in Y_i for the range of magnetic fields found in the far SOL along the ITER divertor, including the highest charge states of neon and oxygen, Be⁴⁺, ⁴He²⁺, and deuterium. ³He²⁺ also experiences a peak in Y_i for magnetic field strengths found within the outer SOL. The results indicate that far-field RF sheaths on divertor surfaces that are not magnetically mapped to the antenna in a device such as ITER would still produce a significant enhancement in sputtering. On the other hand, despite the local minimum in beryllium sputtering yield at resonance (figure 12), off-resonant B-field strengths will play a larger role for beryllium sputtering, especially for the low-field beryllium wall side. This is made apparent by the peak in Y_i within the low-field blue shaded region. Only ³He²⁺ and deuterium have a peak in Y_i for the high-field side beryllium first walls. In short, the maxima for Y_i of both beryllium and tungsten occur in magnetic field ranges where each material will be located in ITER. This demonstrates that the full range of the resonant phenomenon will be relevant for ITER sputtering. Finally, beryllium erosion could become a concern in future tokamaks which will experience longer plasma burn pulses. As ITER is expected to have pulses $\unicode{x2A7E}$400 s, and future tokamaks look to achieve ignition and steady-state, erosion of the beryllium first wall may require further attention and possibly compared against high-Z alternatives producing lower sputtering.

Ultimately, the combination of plasma density, ion temperature, magnetic field strength, and wall material serve to demonstrate locations in a real-world tokamak where enhanced sputtering connected to far-field RF sheaths may occur. The choices for these parameters are reasonable values for the ITER far SOL. Nevertheless, the reader should note our choices of sheath voltage and magnetic field inclination, which, as we have shown, modify the resonance profile. While these values do not affect the device locations that will experience enhanced sputtering, further analysis would be required to quantify the effect of these parameters on sputtering yields. Furthermore, net erosion rates are not addressed by this analysis and await coupling to SOL modelling codes that can refine parameters such as plasma density and temperature.