Time-Domain Modeling of RF Antennas and Plasma-Surface Interactions

Recent advances in finite-difference time-domain (FDTD) modeling techniques allow plasma-surface interactions such as sheath formation and sputtering to be modeled concurrently with the physics of antenna near- and far-field behavior and ICRF power flow. Although typical sheath length scales (micrometers) are much smaller than the wavelengths of fast (tens of cm) and slow (millimeter) waves excited by the antenna, sheath behavior near plasma-facing antenna components can be represented by a sub-grid kinetic sheath boundary condition, from which RF-rectified sheath potential variation over the surface is computed as a function of current flow and local plasma parameters near the wall. These local time-varying sheath potentials can then be used, in tandem with particle-in-cell (PIC) models of the edge plasma, to study sputtering effects. Particle strike energies at the wall can be computed more accurately, consistent with their passage through the known potential of the sheath, such that correspondingly increased accuracy of sputtering yields and heat/particle fluxes to antenna surfaces is obtained. The new simulation capabilities enable time-domain modeling of plasma-surface interactions and ICRF physics in realistic experimental configurations at unprecedented spatial resolution. We will present results/animations from high-performance (10k-100k core) FDTD/PIC simulations of Alcator C-Mod antenna operation.


Introduction
Physical processes occurring in the edge and scrape-off layer (SOL) of magnetically confined fusion plasma can dramatically reduce the effectiveness of the RF sources used for heating or current drive within the core plasma. Sheath formation on plasma-facing antenna components, for instance, may induce sputtering of high-Z ions from these surfaces; such impurity components induce radiative cooling as transport processes carry them toward the reactor core, and may quench the fusion reaction altogether. In addition, parametric decay instabilities or other resonant phenomena may be excited in the plasma edge, reducing RF antenna efficiency as injected RF power is partially diverted away from the core into undesired edge oscillations. In these interactions, the plasma edge density plays a significant role. Higher densities near plasma-facing antenna components are associated with higher material sputtering rates and sheath potential drops, and although reduced densities and sheath potential values can be achieved by widening the gap between the antenna and the plasma's last closed flux surface, such a change may also open additional, undesired channels into which RF wave power may flow.
In this paper, numerical techniques for modeling parasitic wave excitation and the physics of plasma-material interactions in the tokamak edge and scrape-off layer will be described. We will present a self-consistent finitedifference time-domain (FDTD) model describing the concurrent evolution of RF wave fields and sheath potentials, and will describe its implementation in the VSim FDTD code [1]. Half-torus simulations of ICRF heating in C-Mod, which make use of realistic experimental configurations and plasma profiles from experimental data, will be described; these simulations run for hundreds of ICRF periods with millimeter-scale spatial resolution. We will also discuss a process for modeling sputtering impurity generation in these discharges using VSim's particle-in-cell (PIC) capabilities. These complex numerical simulations are possible only due to the increasing capabilities afforded by massively parallel supercomputing platforms, and the computations we describe here exercise these capabilities at large scales (tens of thousands of processor cores).

FDTD methods for EM waves in plasma
VSim uses electromagnetic FDTD methods [2] which are modified to account for the presence of plasma. Electric and magnetic field vector components associated with an individual grid cell are not co-located at the cell nodes, but are evaluated (respectively) on edges and faces of the grid cell in the conventional manner. Leap-frog timestaggering of the electric and magnetic fields relative to one another preserves proper time-centering; the path and surface integrals of Maxwell's equations are evaluated respectively around cell edges and over cell faces to center these equations properly in space.
Reference [3] describes how cold-plasma currents in the Maxwell equations can also be included in the FDTD formalism. Components of plasma current associated with an individual grid cell are co-located at cell nodes [permitting evaluation of the cold-plasma current evolution equation (3); see below]. Interpolation methods which carry edge fields to the grid nodes ("edge-to-node") and nodal fields to edges ("node-to-edge") thus become necessary when such sources are used. A locally implicit time-advance algorithm enables the numerical restrictions imposed by plasma cutoffs (when wave phase velocities diverge) to be circumvented, and preserves the physics of linear plasma dispersion in the presence of a background magnetic field.
Cut-cell techniques [4] are used to model geometrically complex boundaries. Cut edge lengths, associated with portions of individual grid cells intercepted by and lying partially outside conducting boundaries, are used in the integral Maxwell equations. Edge-to-node and node-toedge operations for the corresponding fields must be appropriately modified to reflect the presence of the conducting boundary. At this boundary, the physics of plasma sheath formation can also be modeled.
The extreme length scale disparity between the width of such sheaths (generally a few Debye lengths) and the wavelengths of the ICRF wave phenomena of interest here precludes direct spatial resolution of the sheath. However, salient sheath effects can be represented by a "sub-grid sheath boundary condition" [5], in which local capacitive and resistive circuit elements associated with the presence of the sheath are defined at the nodes of grid cells intersecting the material boundary. This formulation enables a local sheath potential to be defined at the boundary. Reference [6] generalizes the methods of Ref. [3], such that these time-varying sheath potentials at the conducting boundary can be modeled within the locally implicit time advance. Fields at the k-th simulation timestep Δt are evolved to the (k+1)-st timestep by the equations wherein e is the node-centered electric field obtained by applying the non-square edge-to-node matrix operator M to the edge electric field E (with the transposed matrix M T forming the node-to-edge operator), B is the face-centered magnetic field, j α is the plasma current of species α, φ is the local sheath potential, Ω α is the signed species gyrofrequency vector in the background magnetic field, and ω pα is the species plasma frequency; the speed of light c and vacuum permittivity ε 0 also appear. Sheath potentials are defined only at grid nodes with one or more adjoining edges cut by the conducting boundary; vectors normal to this boundary are interpolated to the neighboring nodes outside the boundary to form the unit normal n associated with each such node. The sheath width Δ s and the sheath dissipation parameter ν s are, most generally, spatially and temporally varying numerical parameters defined on these boundary nodes. Equation (3), which computes the cold-plasma response of the species plasma currents, removes the sheath electric field (subtracting the approximate negative gradient of φ across the thin layer Δ s ) to compute this response only as a function of the electromagnetic fields in the bulk plasma. Equation (4) is a Kirchoff relation for the parallel RC circuit whose behavior approximates the physics of the sheath near the wall. Equations (2) -(4) can be solved implicitly.
Explicit FDTD PIC simulations were used in Ref. [6] to quantify the circuit-like behaviors of the sheath. These simulations model the sheath width with the expression which captures DC Debye shielding and Child-Langmuir effects in the limit of small or large 〈φ〉 (compare with Ref. [5]). In Eq. (6), T e and λ De are the electron temperature and Debye length, q i is the elementary charge, and 〈φ〉 is the DC-averaged sheath potential over the RF wave period. More detailed expressions for the sheath capacitance and resistance, such as those of Ref. [7], can also be incorporated into the model, and future code development efforts to incorporate such models in VSim are envisioned.

Alcator C-Mod Simulation Geometry
Simulations are carried out on a rectangular mesh encompassing half of the Alcator C-Mod torus. Beginning with an initially metal-filled domain; toroidally symmetric boundary wall data is read from an equilibrium EQDSK file and rotated to hollow out the 3D vacuum vessel interior. Limiter boundaries corresponding to antenna surfaces are ignored in this data; the approximate shape of the interior vessel with and without antenna limiters is shown in the solid line curves of Figs. 1a and 1b. The CAD file for C-Mod's field-aligned ICRF antenna is imported into the simulation and properly aligned with the formerly axisymmetric wall. Current driven within the coaxial antenna feeds powers the antenna. Open boundary conditions are used at the rear of the coaxial feeds, while a matched absorbing layer boundary condition is used on the non-metallic portions of the plane bisecting the torus. All other boundaries in the simulation are metallic. Magnetic fields are generated from the EQDSK Grad-Shafranov data and are  1050426022 are mapped onto the simulation domain to construct the plasma profile (see Figure 1b). The entire half-torus simulation domain including the plasma (cut away on the left side to reveal the antenna placement) is shown in Figure 2.

Fast and Slow Waves
In Figure 3 we plot a typical wave field generated by the model described above. In this discharge, adjacent antenna straps (of which there are 4 in total) are driven 180 degrees out of phase from one another, yielding complex wave patterns as fast waves are excited in the edge plasma and propagate toward the plasma core with some degree of focusing. Contours in the figure show oppositely signed values of the vertical electric field. In addition to the fast waves, millimeter-scale slow wave excitations in the plasma edge and SOL are visible in the simulations; the slow wave is evanescent outside of a narrow layer near the plasma edge. One focus of present modeling efforts is on determining the fate of the power which flows from the antenna and is diverted into these slow waves, since quantifying power losses associated with parasitic slow wave excitations as plasma edge parameters are varied will enable more efficient antenna operation scenarios to be developed. Fig. 3. VSim modeling of ICRF antenna operation in Alcator C-Mod; contours of vertical electric field (dark/light contours have equal magnitudes and opposite signs) are shown. Fast waves induced by the antenna form periodic wave structures of long wavelength (tens of cm) near the tokamak core; some power is lost in the SOL as millimeter-scale slow waves are excited in the thin layer between the density cutoff and the lower-hybrid resonance. This simulation used 60,000 cores on the Titan Cray XK7, simulating 250 RF antenna periods in 12 wallclock hours.

Sputtering and Impurity Modeling
To model sputtering and impurity physics, VSim uses material composition data for the various antenna surfaces in tandem with its PIC modeling capabilities. Sputtering yields vary depending on the material surface composition, so the model must represent the different plasma-facing surfaces comprising the antenna and vacuum vessel as distinct material objects with correct physical properties. Protective tiling on C-Mod's interior vessel surfaces is made of TZM (a molybdenum alloy with traces of titanium, zirconium, and carbon). Faraday shields (the many nearly-horizontal rods in Figure 4) are also made of TZM. The antenna straps (large metallic strips beneath and perpendicular to the Faraday shields in Figure 4) and limiter structures (grey notched and solid EPJ Web of Conferences 157, 03021 (2017) DOI: 10.1051/epjconf/201715703021 22 Topical Conference on Radio-Frequency Power in Plasmas structures surrounding the straps) are high-purity oxygenfree copper plated respectively over electroless nickel flash on Inconel 625 and 316 stainless steel. VSim imports these different components from within the antenna CAD file as separate physical objects, and associates these with relevant material composition data to model sputtering effects as particles impinge on these distinct surfaces. The particles which induce such sputtering events are treated in the simulations as test particles; they respond to the electromagnetic fields of the background and RF waves but do not contribute to these fields. Particles are probabilistically loaded into the domain such that their density profiles match the experimental edge and SOL profiles. Temperature profiles for these particle distributions are also obtained from C-Mod data. As the simulation proceeds, particles which pass through the sheath to strike the material wall have their strike energy incremented in accordance with the local sheath potential at the strikepoint, enabling the sputtering computations (as well as local measurements of particle and heat flux) to accurately take into account the energetics of sheath transit. We continue to develop the model's sputtering capabilities; current areas of focus include improving ionization models in the SOL for the sputtered neutral wall atoms, as well as improving surface models for the material alloys comprising the wall and antenna surfaces. Figure 4 shows initial data from the sputtering modeling; here, the ICRF antenna was run for many RF periods to determine the locations on the antenna where ion strikepoints were most likely to occur. The possibility that such time-integrated data can be compared with wear/erosion patterns on existing antenna structures in experiments is of interest. In addition, we hope to apply the model to more generally explore sputtered impurity production, impurity transport, and the interaction of RF waves with impurity populations.

Ongoing and Future Modeling Efforts
We have begun efforts to generalize the modeling techniques used here so that they can be applied to other experiments; currently our efforts focus on ICRF modeling for ITER and EAST.
The non-trivial computational resource requirements for this modeling continue to pose challenges; the development of algorithms which make use of modern GPU capabilities to increase speed and efficiency is being explored. Efforts to develop improved sheath boundary conditions, diagnostics, and impurity generation models are also ongoing. In addition to the continuing study of parasitic power losses due to slow wave excitation, we have also begun to explore the physics of density pump-out arising from ponderomotive force effects [8]. The inclusion of yet more detailed sputtering models via coupling to F-TRIDYN [9], a Monte-Carlo Binary-Collision-Approximation code which simulates the physics of plasma-wall interactions at the atomic level, is a topic of interest for integrated modeling efforts. Other potential integrated modeling opportunities include couplings to impurity transport, edge turbulence, or extended MHD codes.