Using the Stix finite element RF code to investigate operation optimization of the ICRF antenna on Alcator C-Mod

As the Ion Cyclotron Radio Frequency range (ICRF) heating becomes more favorable in fusion devices, the urgency of predicting and mitigating impurity generation that arises from it becomes more pressing. In the ICRF regime, rectified Radio Frequency (RF) sheaths are known to form at antenna and material edges that influence negative effects like sputtering and a decrease in heating efficiency. Methods to mitigate the formation of these RF sheaths through RF image currents cancellation have been experimentally studied. A power-phasing scan done on Alcator C-Mod in which the amount of power on the two inner straps (P in) versus the total 4 straps (P tot) was varied showed a minimization of enhanced potentials between Pin/Ptot∼ 0.7–0.9 while impurities were minimized for Pin/Ptot∼ 0.5–0.8. New capabilities in the realm of representing the RF sheath numerically now allow for these experiments to be simulated. Given the size of the sheath relative to the scale of the device, it can be approximated as a Boundary Condition (BC). A new parallelized cold-plasma wave equation solver called Stix implements a non-linear sheath impedance model BC formulated by Myra et al (2015 Phys. Plasmas 22 062507) through the method of finite elements using the MFEM library [http://mfem.org]. It is seen that Stix shows qualitative agreement with the measured C-Mod enhanced potentials.


Introduction
Using the Ion Cyclotron Frequency Range (ICRF) for Radio Frequency (RF) heating is a robust method applied to various fusion devices around the world. Although ICRF heating is * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. efficient, experiments show that in this regime there is a formation of a thin micro-scale RF sheath at the material and antenna boundaries that influences adverse effects such as sputtering, edge power dissipation, and hot-spots [1][2][3]. With future tokamak devices relying on ICRF, like SPARC and ITER, it is becoming increasingly important to numerically model these detrimental RF sheaths through an integrated model approach.
RF sheaths are typically classified into two categories: nearfield sheaths, which form on the antenna and limiter surfaces, and far-field sheaths, which form farther away from the antenna source [4]. The focus of this paper will be on the near-field sheaths that arise when some E ∥ is inadvertently launched by the antenna causing a SW to interact with the Plasma Facing Components (PFCs). These near-field sheaths are highly dependent on the antenna geometry, power, and plasma density near the launcher [4].
There have been numerous studies aimed at mitigating near-field sheaths and the associated impurity generation. These have lead to discoveries such as using field-aligned Faraday screens (FS) [5], low-Z coating on FS [6], using electrically insulating enclosures on antennas [7], using a k ∥ < k 0 = c/ω to reduce propagating coaxial modes and surface waves [8,9], and lastly RF image current cancellation [10][11][12]. Due to the experimental successes on both ASDEX-Upgrade (AUG) and Alcator C-Mod, the current focus is using RF image current cancellation to optimize the antenna geometry and current strap phasing for future devices like ITER and DEMO [13,14]. The main idea behind this method is that by varying the fraction of current in each strap of the antenna, there becomes an optimum fraction that induces enough image current that cancels out the E ∥ from the parasitic SW thereby eliminating the near-field sheath.
In the pursuit of decreasing tungsten (W) impurities during ICRF heating, the AUG team investigated methods of optimizing their two-strap antenna to reduce the E ∥ on the antenna limiters that included reducing RF image currents [3]. This discovery lead to a new design of a three-strap antenna that did see a reduction of W impurity generation [13]. Using this antenna on AUG a straps phasing of 0/π/0, a power-phasing study was done in which the fraction of power of the central strap relative to the outer straps was varied to see if there is an ideal operating scenario [10]. On both sides of the RF limiters, the RF potential (V RF ), the DC current (I DC ), and the tungsten sputtering yield (Y W ) were measured in several poloidal locations [10,13]. It was found that there were variations in the local minima of V RF in power-balance and phasing space on both sides of the limiters [10]. The overall trend indicated a minimum of both V RF and Y W between P in /P tot ∼ 0.55-0.7 with the ideal fraction at P in /P tot ∼ 0.67 [10,15].
Following AUG's E ∥ minimization efforts via image current cancellation, experiments were run on Alcator C-Mod's field aligned four-strap ICRF antenna using the same power phasing method with 0/π/0/π current strap phasing for both P ICRF = 1 and 1.5 MW. Although the C-Mod antenna geometry and SOL density are different, the trends of minimization measured in C-Mod were similar to that of AUG. Using Gas-Puff Imaging (GPI) measuring the upper outer corner of the antenna showed that enhanced potentials were at a minimum between P in /P tot ∼ 0.7 − 0.9 [11,12,15]. In the case of the 1 MW scan, in the region of minimization, the potentials were similar to that of having no RF [11]. The main impurity for C-Mod is Molybdenum and was found to be minimized for P in /P tot ∼ 0.5 − 0.8 for both near and far from the antenna [11,12,15] Although the results from the AUG and C-Mod experiments uncover an ideal operating regime for ICRF antennas, there still is a need to simulate ICRF antenna scenarios in the context of an integrated model that includes both transport effects and the resulting impurity generation. RF sheath potentials tend to have large spatial gradients that give rise to electric fields on the PFC surfaces. When combined with the background magnetic field, E×B flows are created leading to what is known as RF convective cells [16,17]. These RF convective cells are known to influence the density and impurities near the antenna [16,17]. This effect has been seen experimentally on various devices like in Tore Supra [18] and JET [19]. Numerical modeling of RF sheaths is typically done by solving for the wave fields using an Electromagnetic (EM) solver that does not include the transport effects of the RF convective cells which has been shown to be a key influence on plasma-material interactions. Similarly, coupling of RF sheath wave codes to impurity generation codes is lacking. In order to have a complete picture of impurity and RF sheath behavior around the ICRF antenna, an integrated model is vital.
This paper focuses on the initial step towards an integrated model approach of ICRF antennas using the C-Mod powerphasing study as the focus. Presented here is the first component: the RF field solve with the resulting rectified potentials using the Stix code. This RF code was chosen because of it is light-weight nature and its future capability to couple to the transport code 'MAPS' [20]. A detailed description of Stix's internal algorithm can be found at Migliore et al [21]. The overview of the paper is as follows: section 2 introduces the setup of the 2D C-Mod simulations using the RF code 'Stix.' Section 3 discusses the results of a 2D power-phasing scan using Stix. Lastly in section 4, the implications of the simulations and future work are summarized.

Simulation setup using Stix
Stix is a finite-element cold plasma RF wave solver built off the library of MFEM that solves for the magnetic field, ⃗ H, of the launched wave given by equation (1), The plasma dielectric,ε, is given as whereĪ is the identity matrix andb is the unit vector of the background magnetic field defined asb = ⃗ B 0 /| ⃗ B 0 |. Using Stix notation [21], the coefficients of S, P, and D are defined as Within these Stix coefficients, an artificial collisionality was incorporated through an effective mass on the electrons given as m eff = m e (1 + iν e /ω). This is a common method used in EM solvers to add damping in addition to emulating the thermal effects that resolve the singularities of the cold plasma dielectric. These singularities consist of wave resonances and for Stix wave cutoffs due to the inversion of the dielectric tensor.
Stix also solves for the RF sheath through a non-linear finite impedance Boundary Condition (BC) formulated by Myra et al [22,23] given by equation (4).
This sheath BC is expressed in terms of the tangential electric field on the wall, which is taken to be the tangential gradient of the RF sheath potential, V sh , oscillating at the wave frequency. Due to the sheath width being much smaller than the wavelength of the wave, the sheath is treated as a 1dimensional problem and therefore is electrostatic. The RF current passing through the sheath, J n , which includes the electron, ion, and displacement components is characterized by the voltage drop across the sheath, V sh , and a complex sheath impedance z sh [17]. J n in turn can be expressed as the normal component of the wave electric displacement, D n , which yields the third term in equation (4). For use in Stix, the BC is rewritten to be expressed in terms of the ⃗ H field instead of D n . Other RF wave codes typically solve for the electric field, Stix however was chosen to solve for the magnetic field, ⃗ H. This choice was made so that the non-linear sheath BC, equation (4), would be solved simultaneously with the plasma wave equation, equation (1). In the electric field formulation, the entire electric field solution needs to be first solved using a guess for the sheath impedance to update the impedance and therefore the E t boundary term. Using the H field formulation bypasses having to solve the problem in two steps and allows for both equations to be solved together in a block matrix solve leading to faster convergence.
The MFEM library is inherently 3D which results in Stix operating in pseudo-1D and pseudo-2D, where either the y and z or just z direction(s) are represented by phase factor(s) in the form of a plane wave ∼ exp (ik y y + ik z z) and exp (ik z z) respectively [21]. For replicating the C-Mod 4 strap antenna numerically, the simulation domain was taken to be a 2D slice along the 10 degree pitch angle of the background magnetic field in the mid-plane of the field-aligned antenna shown in figure 1. For these simulations, the phase factor that represents the parallel wavenumber, k ∥ , is imposed through the physical dimensions of the four antenna current straps as well as their dipole phasing of 0/π/0/π which resulted in k ∥ = 11 m −1 . The features and dimensions of this 2D mesh were taken from a slice of the 3D CAD of the ICRF antenna that did not include the Faraday screens. The computational domain spanned from R major = 0.4 to 1.1 m and −0.6 to 0.6 m in the tilted quasitoroidal plane. Following the true C-Mod geometry, the front of the RF limiters were placed at R major = 0.913 m and the front of the current straps were placed at R major = 0.935 m. The current straps were simplified to be boxes of current with dimensions of 70 mm in the quasi-toroidal direction and 10 mm in the radial direction in which each strap's total current was set using Dirichlet BC on tangential component of the magnetic field, H t . The poloidal current was assumed to be a cosine squared profile that peaks at the center of the strap near the mid-plane. Given that the poloidal variation is not represented in the Stix 2D plane, the amplitude of the current was taken to be at the peak of the cosine squared poloidal profile that corresponds to the mid-plane of the antenna.
The radial electron density profile was found using a pedestal-like sech profile along the major radial direction adapted from Kohno et al (2015) given by equation (5) where j represents the particle species of deuerium ions and electrons [24]. To get as close as possible to experimental conditions, for this simulation the density profile used two pedestal-like profiles in which λ n,j , ν j , n max,j , and n min,j in equation (5) were parameterized to fit previously measured L-mode density scans from C-Mod [25]. This experimental density data was only measured to the RF limiter which is located at 0.913 m and therefore for this simulation the profile needed to be extrapolated to the vacuum vessel wall using a ∼2.4 mm e-folding depth. The parameterized electron density profile is shown in figure 2. In this profile the lower hybrid (LH) resonance occurs just behind the RF limiter at 0.915 m, the slow wave cutoff (P = 0) occurs just in front of the current straps at 0.934 m, and a density below the plasma frequency occurs at and behind the current straps.
These simulations did not include the high-field side in the domain and therefore to account for no wave reflections, an artificial collisional profile was used through an effective mass of the electrons. This profile was chosen as an exponential profile, similar to damping profile given in [26], which increased in strength as the wave propagates towards the core represented in equation (6) .
Here C 1 represents the strength of absorption, λ 1 represents the e-folding decay length scale, and R min represents the beginning of the domain on the high field side at 0.4 m. Using values λ 1 = 0.035 m and C 1 ∼ 10 2 ω, gave an absorption profile that showed no reflected waves verified through fast Fourier transforms of the resulting electric fields.
In addition to damping out the launched wave, given the density profile used there were two problematic regions that needed to be resolved using artificial collisionality: the LH resonance at S = 0 and the slow wave cutoff at P = 0. Usually cutoffs are not an issue for EM RF solvers, however because Stix solves for the magnetic field, the dielectric tensor,ε, needs to be inverted therefore making P = 0 a singularity. For both P = 0 and S = 0 layers, the artificial collisional profile used was a thin Gaussian represented by equation (7).
Similar to equation (6), C 2 is the strength of absorption, λ 2 is the squared e-folding decay length taken to be 10 −4 m 2 , and R loc is the location of the resonance/cutoff. In contrast to equations (6) and (7) C 2 was on the order of ∼ ω due to the fact that the solver needed only a small amount of finite collisionality to resolve the singularity rather than be strong enough to damp out the wave entirely. For the P = 0 layer, the term with the most influence for these plasma parameters is the electron contribution and therefore modifying the artificial collisional frequency of the electrons, ν e , using equation (7) allowed the cutoff to be resolved. Similarly, the LH resonance can be resolved by adding some finite amount of collisionality to the ion term shown in equation (7).
The last plasma profile needed was the background magnetic field which was taken from an EFIT reconstruction of the magnetic equilibrium of the experimental discharges, stored in an EQDSK file format [27]. This file only has the flux function calculated for the poloidal (R, Z) plane and therefore a transformation into the Stix computational reference frame needed to be done. This was done by transforming the magnetic field values from the poloidal cross section to a 3D torus then finding the values that corresponded to the titled top-down view of the four-strap antenna.
Given the sharp density gradients, the existence of the LH resonance and SW cutoff, and the proximity of these layers to one another, this problem poses difficulty to numerically resolve. Shown in figure 3 is the evolution of the uniform global refinement of the mesh with both linear and second order polynomials on the real parallel electric field using the Stix code. With increasing uniform refinement and 2nd order  Pseudo-color plot of the real parallel electric field for P in /Ptot = 0.5 for various uniform refinement levels (RS) and order of finite element polynomial. The RS value denotes the number of times the mesh has been refined uniformly with 0 as no refinement. This region is highlighting the slow wave propagation between the front of the current strap boxes and the RF limiter that is responsible for this rectification. The noise is due to the sharp density gradient in addition to the LH and P = 0 SW cutoff nearby which needed to be resolved with artificial collisionality. polynomials, the solution of the SW propagating in this thin region becomes increasingly visible. The first and second most resolved cases showed the H field solution's norm varied by 0.01% from one another deeming the solution to be resolved. Due to Stix being inherently three dimensional, the number of uniform refinements is limited to the amount of available computational memory. Every uniform refinement increases the number of elements by 4 × therefore allowing for only a few uniform refinements. The most intensive simulations were the order 2 polynomial with uniform refinement of 1 with ∼5 × 10 5 elements at ∼1.7 × 10 6 number of unknowns. From these refinement simulation studies, a key takeaway is that in the region in between the RF limiters and the front of the current strap using a realistic density profile, it is necessary to have significant resolution using higher order polynomials along with some non-negligible artificial damping at the S = 0 and P = 0 layers when using a cold plasma model.

Resulting simulated rectified potentials and E ∥
For the Stix's power-phasing simulation scan, the sheath BCs were placed on the faces of both the left and right RF limiters, shown by the orange lines in figure 1(a), bypassing regions in which the magnetic field is tangent to the surface where the sheath BC breaks down [22]. All other domain boundaries for these simulations were taken to be conducting wall in which lossy materials were not considered. Emulating C-Mod's experiment, a pure deuterium plasma was used along with a background magnetic field of 5.4 T on axis read in from the EQDSK. A constant temperature profile of 10 eV was taken which was only used within the sheath BC. Using a 80 MHz launched wave, power-phasing fractions of P in /P tot = 0.05 to 1.0 were scanned using a 0/π/0/π antenna strap phasing. Both 1 and 1.5 MW cases were run in which the total power of the four straps was held constant and the total current of each strap was adjusted to match the given power fraction. The resulting parallel electric field and rectified potentials were measured on both the left and right RF limiters and found to have the largest values on the inner facing components. This region of enhancement aligns with the fact that the C-Mod density profile traps the propagating SW in the region between the front of the straps and the RF limiter. Plotting the real parallel and perpendicular electric fields in figure 4 shows the propagating SW and its intersection with the limiter corresponding to where the rectification was found.
Shown in figure 5 is the resulting DC rectified potentials (V REC ) taken at the location where the largest potential was found. This peak occurred on the left inner facing limiter for both 1 and 1.5 MW cases. Similarly, figure 6 shows the parallel electric field magnitude taken at the same point as in figure 5. In both power-phasing scans, one sees that there is a minimization of E ∥ and V REC for P in /P tot ∼ 0.8 − 0.95. The only difference in trends between V REC and E ∥ is that there is less of an increase towards the highest fraction in the V REC plot. This distinction is due to the fact that the RF sheath rectification is a non-linear effect which was also observed in TOPICA runs of the JET A2 antenna [15]. Using a linear sheath BC, taking z sh ∼ ∞, would show the rectified potentials follow the same trend as the parallel electric field and would overestimate their values particularly towards the highest P in /P tot fractions.
Both figures 5 and 6 show that Stix's 2D reference frame along the background magnetic field in the middle of the antenna does not produce a large amount of rectification given that the Bohm sheath is at 31.86 V. This suggests that there are 3D geometrical effects not accounted for. The importance  of these 3D effects is supported by the experimental measurements and TOPICA simulations of V RF and the tungsten sputtering yield (Y w ) found to significantly change in strength and shape in the poloidal direction on the limiters in the AUG 3-strap antenna [10,15]. Comparing the values of the measured V RF on AUG in the region closest to Stix's reference frame shows that V REC are at the lowest values agreeing with the low voltages found for Stix's scans. In addition to the 2D reference frame, the amplitude of V REC is strongly dependent on the strength of the SW. While not as realistic to the true density in the experiment, it was seen in Stix simulations using a density profile that pushed the SW cutoff behind the current straps resulted in ∼20 V higher values of V REC . Even though changing the density profile resulted in a different amplitude of V REC , the minimization stayed at the same power fraction range. This effect was also seen in similar TOPICA and SSWHICH simulations of the JET A2 antenna [15].
Comparing to the C-Mod measured enhanced potentials given in [12], Stix's V REC show qualitative agreement in the trend that minimization occurs at around P in /P tot ∼ 0.85. It is difficult to compare the GPI measurements of the plasma potential to the V REC found in Stix due to the fact that they are measured in different locations and that the experimental data measured was the E×B velocity taken to be proxy for the DC potentials. The simulated enhanced potentials from Stix are taken from the inner facing left RF limiter near the midplane of the antenna while the experimental data is taken in the upper outer corner of the antenna box. It has been seen both experimentally and in EM simulations that there is a larger enhancement of electric fields and rectified potentials around the corners of the antenna box than at the mid-plane of the antenna that would explain the difference in rectification strength of Stix versus experiment [10,28]. This difference in where V REC was measured between simulation and experiment would also explain the shift of the minimum towards a higher power fraction seen in Stix. This effect of the poloidal variation in the minimum of ⟨V RF ⟩ was seen in TOPICA and Raplicasol simulations that showed a shift towards the highest fractions in the same regions as the Stix 2D reference frame [10]. More generally, due to the complicated nature of this region in the plasma it is difficult to get a quantitative analysis of V REC for comparing computational and experimental data. There are two notable limiting factors: (1) one cannot truly know the density in front of antenna straps since it is inaccessible to density diagnostics and (2) there are transport effects like the E×B flow that are neglected. Transport effects play a key role in RF sheath rectification as discussed in the introduction [16]. These effects can be seen experimentally in the three-strap AUG antenna which measured density fluctuations in various locations along the front of the antenna that would affect the resulting enhanced potentials [29].
Lastly, it should be noted that previous experiments on the field-aligned C-Mod antenna using monopole (0/0/0/0) phasing showed a significant increase in enhanced potentials [30]. Due to all the current strap's phasing being in the same direction, the image current cancellation is drastically reduced allowing for much more sheath rectification. Using the same plasma conditions and background magnetic field, running Stix's C-Mod configuration with monopole phasing for P in /P tot = 0.5 to 1.0 confirmed what is seen experimentally. A comparison of the peak sheath rectification for the monopole and dipole phasing is shown in figure 7. Figure 7(a)'s potentials are taken at the largest point of rectification on the right RF limiter's inner side in contrast to figure 7(b)'s which is taken on the left RF limiter's inner side for both phasing regimes. On both sides of the antenna's RF limiters, the rectification is enhanced with monopole phasing with no minimization. Figure 7 also shows an asymmetry of the strength of the sheath potential, an effect that is seen in both monopole and dipole phasing Stix simulation scans but which is exacerbated with monopole phasing. This effect comes from the poloidal magnetic field variation and is also seen experimentally as shown in [10].

Summary and conclusions
This paper introduces the first step of an integrated model approach for numerically simulating RF sheath physics on ICRF antennas using the finite-element EM solver Stix. The focus of this study was to qualitatively reproduce the minimization trend found on the Alcator C-Mod four strap ICRF antenna power-phasing study. It was found experimentally that the enhanced potentials were minimized between P in /P tot ∼ 0.7 − 0.9 while the Molybdenum (Mo) impurity flux was minimized between P in /P tot ∼ 0.5 − 0.8 for both close and far from the antenna [11,12]. Using the reference frame of a 2D field-aligned toroidal slice along the mid-plane of the C-Mod antenna, Stix found that V REC and E ∥ minimized between P in /P tot ∼ 0.8 − 0.95 giving qualitative agreement with the experiment.
The results found from Stix and other EM RF simulations of minimizing E ∥ show that this is a complex problem. It was seen that in the AUG experiment and the consequent TOPICA simulations that the minima of the enhanced potential varied with poloidal location on the RF limiters [10]. This suggests there will be regions of rectification that cannot be fully cancelled out with image current. Additionally, TOPICA simulations of the JET A2 and ITER antenna show an interesting effect of how different the minimization trends are for the different antenna geometries [15]. The JET A2 antenna minimization was found to occur at P in /P tot ∼ 0.65 − 0.8 while the minimization on the ITER antenna was shifted to higher fractions P in /P tot ∼ 0.85 − 1.0 [15]. These results emphasize that antenna geometry is expected to play a role in the minimization power ratio as that ratio depends on the relative spacing of the straps and distance to limiter surfaces and why it is important to continue simulating the enhanced potentials on various antennas in a more complete approach.
While minimizing E ∥ on the antenna is important for decreasing the enhanced potentials, the mechanism behind the behavior of the associated impurity generation measured in the various experiments is still not fully understood. This effect can be highlighted by comparing the impurity results found on AUG and C-Mod that show the importance of geometry. In C-Mod, Mo sources were more prominent farther away from the antenna than in AUG because the RF limiter is behind the plasma limiter. In contrast, the AUG RF limiters are the plasma limiters and therefore the tungsten is more local as a source [15]. Not only does the antenna geometry have influence, as mentioned before, transport effects are important in this region and are completely missing with EM RF simulations. The impact of transport can be seen by comparing the minimization region of the Mo impurities versus the enhanced potentials measured on C-Mod in which the impurities were minimized for a lower power phasing fraction and for a broader region [11,12].
The ultimate aim of simulating this power-phasing study of the C-Mod antenna using Stix is to lay down the framework of an integrated model. The two physical processes that would be key to providing a more comprehensive model of realistic plasma conditions around the antenna are a transport model using the code MAPS [20] and an impurity flux calculation using the code RustBCA [31]. Currently, there are efforts underway to couple Stix and RustBCA so that initial impurity fluxes may be obtained from the calculated rectified potentials. These efforts will be reported on at a later date.
phasing study and the implications studies like these have on resolving near-field sheaths. This work was supported by the U.S. Department of Energy Scientific Discovery through Advanced Computing Initiative, Contract Number DE-SC0018090.