Recent progress in modeling ICRF-edge plasma interactions with application to ASDEX Upgrade

This paper summarizes recent progress in modeling the interaction between ion cyclotron range of frequency (ICRF) waves and edge plasma with application to ASDEX Upgrade. The basic theories, the development of ICRF and edge plasma codes, the integrated modeling methods and some key results are reviewed. In particular, the following physical aspects are discussed: (1) ICRF power coupling; (2) slow wave propagation; (3) ICRF-rectified sheath; (4) ICRF-induced convection; (5) ICRF-edge turbulence interaction. Moreover, comprehensive integrated modeling strategies by including all necessary codes in one package and solving multiple physical issues self-consistently are discussed.


Introduction
Radio-frequency (RF) heating with waves in the ion cyclotron range of frequencies (ICRF) is a well-established heating method on present-day magnetic confinement fusion devices, including tokamaks and stellarators. It is also considered as * Author to whom any correspondence should be addressed. a Dr A. Messiaen  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. one of the main auxiliary heating techniques in ITER, where the ICRF heating system is expected to provide 20 MW (baseline) to 40 MW (upgraded scenario) heating power. Understanding the plasma-wave interaction is crucial for the success of this auxiliary heating system. However, the physics involved are usually quite complex, with inter-connections between the various mechanisms in place. Thus, modeling by coupling multiple RF and plasma codes is necessary to describe the physics, to solve the problems and to shed light on future development of the heating system.
Previous efforts have been devoted to understanding and improving ICRF heating. Reviews of the ICRF physics and technology in support of ITER can be found in [1][2][3][4][5][6][7][8]. Progresses in core ICRF modeling and cross benchmarking of RF codes are discussed in [9][10][11]. The development of RF sheath theories and modeling is reviewed in [12][13][14][15]. Some advances of integrated RF modeling are discussed in [16][17][18]. More recently, the RF sheath physics is reviewed in [19] and the ICRF-induced wave-SOL interaction is reviewed in [20]. In this paper, we will focus on discussing the progresses made in recent years in modeling of ICRF-edge plasma interactions. The modeling strategies and some of key results will be discussed.
Modeling of ICRF waves and their interactions with the plasma and material surfaces begins at the antenna where the coupling of power into the plasma is the primary goal. This is also the region where the strongest RF electric fields are normally expected. These RF fields can interact with the plasma and the antenna hardware through the formation of RF sheaths which can potentially lead to both parasitic power dissipation in the scrape-off layer (SOL) and the release of impurities into the plasma.
Since 2012, a self-consistent RF sheath model has been developed by coupling the RF and direct current (DC) parts with non-linear RF and DC sheath boundary conditions at both ends of open magnetic field lines [21]. Based on this model, a multi-2D self-consistent sheath code SSWICH [21,22] was developed to calculate the RF voltage, the DC plasma biasing voltage and the sheath capacitance self-consistently. The recent development of the SSWICH code allows it to handle contributions both from the slow wave (SW) and the fast wave (FW) [23]. More recently, the study of RF sheath theory and models have been significantly developed by parameterizing the RF sheath impedance, the sheath-rectified DC potential and net DC current flow through the sheath [24][25][26] and using oblique magnetic field lines at the sheath entrance [27]. By incorporating this sheath impedance model with the 2D selfconsistent finite element code rfSOL [28,29], the interaction between waves in plasmas and sheaths was studied [29]. In addition, by coupling a 1D edge plasma and sheath code with the plasma-material interaction code Fractal-Tridyn [30], the enhanced sputtering yield from the Faraday screen and the plasma facing components of an antenna can be calculated [31,32].
ICRF waves interact with the edge and SOL plasma in other ways. RF sheath-induced plasma convection and turbulence form a system in which the plasma profiles and RF wave antenna-plasma coupling and wave propagation can interact. Specifically, ICRF wave fields drive plasma drifts and change the edge plasma density while the edge plasma rectifies the RF near-antenna electric field in return. A 1D model has been developed to self-consistently calculate the diamagnetic drifts, E × B drifts and ponderomotive drifts near an ICRF antenna by solving a wave equation, an equation of motion and a continuity equation both on a fast time scale (RF fields) and a slow time scale (macroscopic density and flows) [33]. However, the 3D nature of ICRF-driven drifts necessitates 3D simulations. By using the experimental measured potential in the 3D edge plasma fluid and neutral particle code EMC3-EIRENE [34,35], or by iteratively running the EMC3-EIRENE, RAPLICASOL [36,37] and SSWICH [22] codes in a selfconsistent way [38], the ICRF induced E × B drifts as well as the diamagnetic drifts were successfully calculated. The 3D ponderomotive drifts are planned to be included in these simulations in a way similar to the E × B drifts. Similar work to assess the role of ponderomotive effects in various RF frequency regimes is underway using the VSim code, also discussed below [39]. RF ponderomotive effects have also recently been modeled and compared with measurements on a linear device [40].
The RF waves influence edge plasmas while the edge plasma can also influence the RF waves and power coupling. The amount of power coupled from the antenna to the plasma critically depends on the width of the FW evanescent layer in front of the antenna, which is determined by the plasma density profile in the SOL. The use of the local gas injection to influence the density profiles in the SOL in order to reduce the evanescence layer and, thus, to maximize the ICRF power coupling further motivates a quantitative characterization of the 3D SOL density profiles and antenna coupling resistances. Integrated modeling by coupling the 3D edge plasma code EMC3-EIRENE with the antenna codes such as FELICE [41], ANTITER [42] and RAPLICASOL [36,37] was performed for ASDEX Upgrade (AUG) [43,44], JET [45] and ITER [46]. To study the influence of resonant magnetic perturbation (RMP) fields on SOL density and ICRF power coupling [47,48], the PARVMEC code [49,50] and BMW [51] codes have been used. These compute the 3D MHD equilibrium in RMP discharges including the SOL magnetic field, which is then used in EMC3-EIRENE for the reconstruction of the resulting 3D density profiles. Such 3D density is input in RAPLICASOL for the computation of the antenna S-matrices. Besides, the SOL turbulence can scatter the RF waves and influence the heating efficiency (this effect is in general more significant for lower hybrid and electron cyclotron range of frequency waves owing to their small wavelength relative to that of turbulence), while the RF generated convective flow in the SOL can in turn influence turbulence, especially the filaments. These physics were studied [52,53] by coupling the MHD code JOREK [54,55] or turbulence code BOUT++ [56] with the 3D antenna code RAPLICASOL. An additional mechanism for SOL power dissipation due to ICRF-turbulence interaction is that the incident long-wavelength FW resonantly excites a short-wavelength slow mode at the surface of filaments. More details can be found in sections 2.5 and 4.5.
Progress has also been made in developing antenna and ICRF heating codes. For instance, a 3D finite difference time domain (FDTD) program VSim [57] has been developed to calculate RF wave fields near the antenna, and was successful applied to simulations for Alcator C-Mod and ITER [58,59]. The FW2D code [60] is developed to solve the cold plasma full wave equations to calculate the SOL losses of HHFW in NSTX/NSTX-U. A reduced 2D axisymmetric finite element model (FEM) resolving arbitrary tokamak geometry was developed to study helicon power loss to the SOL [61]. These developed codes, though with reduced physics, can correctly address some specific issues in an efficient way.
To include multiple RF physics from the SOL to the core, such as slow wave propagation and lower hybrid resonance, RF sheath and impurity sputtering, ICRF power coupling and edge power loss, RF induced drifts (ponderomotive drifts, E × B drifts, diamagnetic drifts) and convection, wave scattering by turbulence, wave absorption and mode conversion in the core, current drive, fast ion generation, influence on core plasma transport and turbulence, efforts are required to develop and integrate different codes/models. The important physics in the plasma edge that need to be considered are also shown in figure 1. In this paper, we will focus on discussing the recent progress of integrated modeling to solve several important issues. The eventual goal of the work reported on here is to provide the code building blocks or modules that include the physics of the RF-plasma interactions in the SOL can be used in future large-scale integrated modeling efforts. Discussion of large-scale integration itself is beyond the scope of the present paper.
The rest of paper is organized as follows. In section 2, the theoretical basis of several RF physics topics is discussed. In section 3, the antenna codes, sheath codes and edge plasma transport/turbulence codes are discussed. These codes are either involved in solving various RF issues or related to integrated modeling in the following sections. In section 4, modeling of ICRF-edge plasma interactions is discussed, with emphasis on ICRF power coupling, ICRF induced sheath and convection, ICRF-edge turbulence interaction and coupling of edge-core simulations. Finally, discussions of integrated modeling strategies are given in section 5 and conclusions are given in section 6.

ICRF wave modes
At frequencies typical for ICRF, solving of the wave equation in cold magnetized plasma results in the two wave modes-FW and SW [62]: cs ) are the Stix parameters, Ω cs = q s B m s is the Larmor frequency in plasma with magnetic field B for species s of plasma particles with mass m s , charge q s and density n s , and ω ps = q 2 s n s m s ε 0 is the plasma frequency for particle type s and ε 0 is the permittivity of free space. The wave modes here refer to plane waves oscillating as exp(ik 0 (n ⊥ x + n z)) in the Stix frame for homogeneous plasmas, where k 0 = ω/c while n ⊥ and n are the components of a refractive index.
The lower hybrid resonance (LHR) for SW is defined by the formula and is relevant for the wave equation solution corresponding to the SW. At densities of LHR (10 16 -10 17 m −3 ) the SW goes from propagation to evanescence.

ICRF power coupling
ICRF relies on the FW to transport the wave energy from the antenna at the plasma boundary to the plasma core. However, in the plasma edge, especially in the far SOL right in front of the antenna, there often exists an evanescent layer, because the plasma density is below the FW cut-off density (usually of the order of 10 18 m −3 for typical frequencies and antenna spectrum in nowadays tokamak). For a given anti-node voltage of the transmission line (V max ) and the transmission line characteristic impedance (Z c ), the coupled ICRF power depends exponentially on the evanescent distance (d evan ) from the antenna to the cut-off density, and can be expressed by where R c is the coupling resistance, α is the tunneling factor and k is the parallel wave vector [63]. The above empirical formula of R c is derived analytically under some assumptions [64] while the α value is often obtained from an experimental database. A more strict way to calculate R c , which is often adopted in simulations, is the following [35]: in which Z 0 is the input impedance, SWR s is the standing wave ratio and Γ s is the reflection coefficient for port s. FS i is the feeding scheme and S s,i is the scattering matrix. The above calculations of R c is also used in experiments, where the reflection coefficients can be directly measured by the directional couplers in the transmission line. RF simulations have recently tried to reproduce the whole scattering matrix (or equivalently the input impedance matrix) of multi-port ICRF wave launchers, including the self-inductance of current straps and the inter-strap coupling (mutual inductance) [37]. In principle, the scattering matrix summarizes all the electrical properties of a given antenna face to a given plasma, i.e. it describes the RF system more completely than a coupling resistance. Reproducing the smallest elements of the scattering matrix puts strong constraints on the RF simulations. More details on modeling of the scattering matrix can be found in [37] and measurements on JET can be found in [65].
By increasing the plasma density in front of the antenna, the evanescent distance can be made smaller and the ICRF power coupling can be increased. The edge plasma density can be effectively increased by puffing the fueling gas locally in the main chamber of the machine. Recently, many experimental and numerical studies have been made of the maximization of ICRF power coupling with local gas puffing in many devices, such as AUG [43,44,[66][67][68], JET [45,69,70], DIII-D [71,72], EAST [73,74], WEST [75], ITER [46] and DEMO [76]. It is found that compared to divertor gas puffing, midplane gas puffing close to the antenna can increase the coupling resistance by ∼120% and top gas puffing by ∼20%-40%.
No clear degradation of plasma energy confinement was seen for the applied gas puff rate (on the order of 10 21 -10 22 electrons/s). The SOL plasma density is significantly and locally increased in front of the antenna during midplane gas puffing, and is only moderately and toroidal evenly increased during top gas puffing. Moreover, the impurity sputtering from the wall are found to be decreased during local gas puffing in AUG [68], which is thought in the first place due to the local decrease of edge plasma temperature [77] among many possible effects. In ITER, the ICRF heating system is expected to provide 20 MW (baseline) to 40 MW (upgraded scenario) heating power. However, in ITER d evan can be very large because the antenna is located far from the high density plasma; for example it is ∼0.12 m for a standard lower density profile [46]. Large coupling distances will result in a low ICRF power coupling which is a particular concern in fusion reactors where plasma material issues force a large voltage stand-off. Thus, it is very important to find strategies to reduce the coupling distance such as by tailoring the SOL density profiles with local gas puffing, to maximize ICRF power coupling in ITER.
While maximizing ICRF coupling of power into the plasma is clearly important, increased plasma density at the surfaces of the antenna and nearby hardware can result in enhancement of unwanted plasma material interactions there. RF sheaths, discussed in the next section, are believed to be an important mechanism for these deleterious interactions. Optimization of power coupling while keeping edge interactions, such as parasitic power loss and impurity sputtering at acceptable levels, is one of the challenges requiring interpretive and predictive integrated simulation modeling.

RF sheath
A Debye sheath, whose width is a few Debye lengths, is a thin layer where the plasma interacts with a wall. It is developed because more electrons lost to the wall than ions owing to the much larger mobility of electrons. As a result, it has a large population of ions than electrons. The sheath width is typically much shorter than the other length scales of interest, resolving the sheath is not feasible in codes which include the whole antenna geometry. Thus, in such codes, the sheath must be approximated as a boundary condition (BC). If a large parallel electric field E , often corresponding to the SW, is generated near the wall, it accelerates both the ions and electrons lost to the wall along open magnetic field lines. The more mobile electrons are lost to the wall much faster than the ions. To retain ambipolarity, the plasma develops a large sheathrectified potential to repel most of the electrons. As a result, the sheath is rectified and the potential along these magnetic field lines is enhanced. Understanding the RF sheath is important to solve issues such as RF-related impurity sputtering, parasitic power dissipation and hot spots. Thus, modeling of RF sheaths is required in integrated simulations for a complete description of the effect of ICRF on the discharge. Considering an electrostatic model for the oscillating sheaths, the RF sheath BC for a capacitive, 'immobile ion' sheath can be described by [78]: where the subscripts 't' and 'n' denote tangential and normal (to the wall) components, respectively. α is the incident angle of the magnetic field with respect to the wall. V RF is the instantaneous oscillating voltage across the sheath and V DC is the rectified DC sheath voltage. D n is the electric displacement normal to the wall, δ sh (α) is the sheath width and ε sh is the sheath dielectric permittivity. Recently, a more advanced way to calculate V RF , taking into account particle as well as displacement currents, has been developed by employing the complex sheath impedance parameter z sh [24,27]: where z sh is a function of four dimensionless parameters evaluated on the sheath surface, namely the ion mobility ω/ω pi0 , the ion magnetization Ω i /ω pi0 , the magnetic field direction relative to the surface b n and the RF amplitude eV RF /T e . In this formalism the current density normal to the sheath surface is simply −iωD n . Besides the above dimensionless parameters mentioned, the normalized DC current density through the sheath is further added [24]. For easier implementation in 'global codes', references [24,25] offer simple analytical parametrizations of z sh as a function of the dimensionless parameters, that were validated against 'Debye-scale' codes. In a complement formulation the particle current density normal to the surface is represented as [23]: (13) in which J ⊥B is and J ⊥B es are the ion and electron saturation currents, respectively. The superscript ⊥B represents the situation where the wall is perpendicular to B 0 . V f is the floating potential and V b is the biasing potential expressed as the modified Bessel function, i.e. V b = ln I 0 (V RF ). The above discussed sheath BC and current density equation are often used in sheath codes to calculate the RF sheath-rectified potential, for instance in the SSWICH [22] and rfSOL [28] codes. The philosophy these codes is to split the non-linear electrodynamics in presence of RF-sheaths into RF and DC parts coupled by non-linear sheath electrical characteristics. Both RF and DC sheath BCs are therefore required to solve the coupled model. Equations (11) or (12) apply to the RF part while equation (13) applies to the DC part. The RF/DC splitting is evident in SSWICH but less apparent in rfSOL, where floating DC part is implicitly assumed (i.e. J n = 0). The sheath BCs match the RF plasma current normal to the surface to the current flowing through the sheath. In addition, the RF sheath BCs replace RF quantities evaluated at the metallic walls with RF quantities E t and D n evaluated at the sheath-plasma interface, in which the tangent RF electric field E t is the tangential gradient of the sheath oscillating potential V RF .
Since the RF enhanced sheath-rectified potential can accelerate ions entering the sheath and cause impurity sputtering, large efforts were devoted to minimize this RF sheath rectification. To this end, AUG developed the three-strap antenna, which can significantly decrease impurity sputtering and impurity concentration when operating with optimized phasing and power ratio between straps [79]. Alcator C-Mod used the magnetic field aligned antenna, which also showed a reduction of the impurity contamination and antenna impurity sources [80]. Insulating limiters, which are expected to reduce the induced image current on the limiters and the consequent parallel electric field, have been employed in NSTX/NSTX-U [81], Phaedrus T [82] and Alcator C-Mod [83] to reduce sheath voltages (not completely), but the proper choice of material for extrapolation of that solution to a reactor environment is not apparent.

RF induced drifts
RF-induced DC drifts impact the evolution of the plasma profiles, in particular the density profile, which is important for coupling, sheath interactions, wave propagation and turbulence. ICRF mainly induces two drifts in the edge plasma: the E × B drifts caused by the spatial inhomogeneity of enhanced sheath-rectified potential and the ponderomotive force driven drifts caused by the spatial inhomogeneity of electric field. They are expressed by the following formulas: Here, V DC is the RF rectified sheath-rectified potential and Φ pond is the ponderomotive potential. For electrons and ions, the ponderomotive potential can be written as [33]: in whichα is a factor of order 1 and accounts for the difference in magnitude between cyclotron frequency and plasma frequency for the polarization [33]. Thus, for the ponderomotive drifts, the electrons are dominantly influenced by the parallel electric field while the ions are primary influenced by the perpendicular electric field. Typically, in the ICRF regime the ponderomotive force from FW fields is negligible; however, near the antenna, the SW contribution may be significant.
For the E × B drifts, the electrons and ions feel the same force and they possess the same density convection. These drifts mainly exist on the poloidal cross-sections since the sheath-rectified potential is almost constant along field lines, and results in the so-called convective cells as the RF rectified potential structures are filamentary, and in that respect similar to blobs. For the ponderomotive force driven drifts, they also possesses a non-negligible component parallel to B. Equations (16) and (17) show that this force tends to expel the charged particles from the zones of large RF electric field amplitude. This parallel expulsion likely depletes the density in front of the ICRF antennas. Besides the ions and electrons feel different ponderomotive forces. The ponderomotive transport being non-ambipolar, the electric potential V DC tends to adjust locally in order to equalize the (divergence of ) ion and electron fluxes and thereby preserve the local electro-neutrality of the plasma outside the sheaths. Therefore V DC may not be 'almost constant along field lines' in presence of ponderomotive effects.
Besides these ICRF induced drifts, there are several plasma drifts that originally exist in the SOL, such as the diamagnetic drifts due to the gradients of plasma pressure and the E × B drifts caused by the inherent plasma potential. Since this plasma potential (∼10 V) in the SOL is much smaller than the RF sheath-rectified potential (∼10 2 -10 3 V), the E × B drifts caused by ICRF is much larger than by the plasma itself. In fact, the ICRF induced E × B drifts and ponderomotive drifts in the SOL can be of the same order as the diamagnetic drifts [84]. More generally, the E × B and ponderomotive forces are presented as two independent processes that add-up their effects. The momentum balance could as well be seen as a global problem coupling the E × B, ponderomotive and diamagnetic components, together with parallel pressure gradients. This could possibly be addressed via integrated modeling.

RF wave scattering by blobs/filaments
Within the SOL plasma, poloidally and radially localized density perturbations, called 'blobs' or 'filaments' are ubiquitous [85]. Large filaments (large size and large change in density) occur during ELMs, but the small filaments are always present, even in ELM-free scenarios. RF wave propagating through the SOL can be scattered by filaments, even for ICRF FW which has a much larger wavelength than the filament.
Previously, Ram and Hizanidis derived a general solution for plane wave scattering at a filament within cold plasma theory, valid for all frequencies and wavelengths [86]. Valvis generalized the solution to the case where the filament is not aligned with the confining magnetic field [87], and studied the force the incident wave exerts on the filament [88]. Zhang generalized the result to Gaussian filaments [89]. Myra recognized the possibility that wave-filament interactions can be resonant [90], albeit under the rather unusual condition of a filament with a density lower than that of the surrounding plasma [91]. Tierens argued that for the case of high harmonic fast waves, such as those used in NSTX, the incident fast wave resonantly excites surface waves on the filaments [92], without requiring any unusual filament properties. This is possible because of the excitation of the much shorter slow wave due to BCs on the filament-background interface. The case of non-idealized filaments has been studied in [89,[93][94][95]. One condition under which the FW resonantly excites SW is high harmonic fast wave (HHFW) heating or current drive [96], when is not much greater than 1. Here, r f is the filament radius, K and I are Bessel functions evaluated at k r f m i m e , and ω Ω i 1 for HHFW.
Full 3D numerical treatments of ICRF wave-filament interactions are inherently difficult due to the need to resolve both the slow wave length scale, the density gradient length scale at the filament-background boundary, and possibly very short azimuthal length scales [93]. There are also open questions regarding appropriate BCs with which to 'terminate' the filament in the toroidal direction. A full tokamak geometry, where the filament can be traced all the way to the divertor, is rarely feasible to simulate, and even then, perfect electric conducting BCs at the divertor would not be appropriate [19]. The use of perfectly matched layers [97] to toroidally terminate the filaments, as in [98], requires the surface waves to be toroidally forward [94,99].

ICRF-edge: modeling codes
In this section, various ICRF antenna codes, sheath codes and edge plasma transport/turbulence codes are introduced. These codes are either used in the work discussed in section 4, or are considered as candidate codes for full integrated ICRF modeling discussed in section 5. For each code, the basic equations, physics/mechanisms, recent development/validations, advantages and possible limitations are discussed. The codes discussed here exemplify the present status of modeling, but are not intended to constitute a complete list. For the antenna codes described in section 3.1, the main inputs include density profile, temperature profile (if include hot plasma model), concentration of minority ions, magnetic field, wave frequency, antenna geometry, and the amount of power or voltage given on the port. The outputs usually include magnetic and electric fields, strap currents and scattering matrix. For the sheath codes described in section 3.2, the main inputs include density profile and electric fields. The key output is the sheath potential. For the edge plasma fluid/MHD codes described in section 3.3, the main inputs are the plasma equilibrium, the first wall geometry, the heating power and radiated power, the particle sources, the reference experimental density and temperature profiles. The main outputs are the 3D distribution of the plasma density, plasma temperature, neutral density, Mach number and so on. Comparisons of the main features of different antenna, sheath and edge plasma codes are shown in the tables 1-3.   [100,101] is a combination of two codes: one method-of-moments code that models the electromagnetic fields near the 3D antenna, which is assumed to be in vacuum. This is necessary, since in vacuum the Green function for the electromagnetic wave equation is known, which makes it possible to use a method-of-moments approach. Consequently, a drawback is that TOPICA cannot assess RF sheaths on the antenna structure. The advantage of the methodof-moments approach is that the complicated 3D antenna geometry only requires meshing on the surfaces. Meshing of the antenna volume is, in stark contrast with finite element approaches, not necessary. The degrees of freedom are only the surface currents on the antenna, not the fields in the vacuum volume. The second code is a 1D hot plasma code, FELICE, which represents the fields as functions of the minor radius coordinate 1D and a Fourier expansion in the toroidal and poloidal directions. The coupling of these two codes is the computational bottleneck of TOPICA. It requires the solu-tion of a dense system of equations, whose size is the number of degrees of freedom on the interface (surface) between the antenna in vacuum and the plasma. TOPICA has been thoroughly validated against experiments [102], and continues to be used for the simulation of future ICRF antennas including the ITER antenna [103].

RAPLICASOL.
RAPLICASOL is a finite element code which solves Maxwell's equations in frequency domain in cold plasma. It is implemented within the commercial finite element software COMSOL. Thus, it is very flexible to include complex antenna geometries and build computational grid with different grid resolution at different domains. Previously, it was thoroughly validated against TOPICA for the AUG antennas [37], and more recently also for the ITER antenna [104]. RAPLICASOL does not model the core plasma: the simulation domain is terminated at the plasma side by a layer of fictitious absorbing material, a 'perfectly matched layer' (PML) [97], which mimics the core plasma absorption.
Recent developments of RAPLICASOL include more detailed antenna geometries such as a curved model of the AUG two-strap antenna and the ITER antenna, and the ability to handle 3D density profiles, which have allowed us to model the influence of gas puffing and magnetic perturbations on the power coupling [44,47]. Also, attempts have been made to model the extreme edge plasma where the density is low enough for the slow wave to propagate, which is possible only if collisions are properly included [105].
is a finite element code that solves Maxwell's equations in cold plasma in the edge plasma, around a realistic antenna geometry. PETRA-M can be coupled to TORIC, a hot plasma code that can model the core plasma. A unique feature of PETRA-M is that it can use finite element basis functions of arbitrary order: fourth and fifth order has been used in practice for comparisons with RAPLICASOL, and up to 12th order has been tested. This enables PETRA-M to exploit a trade-off between element order and number of degrees of freedom, enabling it to give accurate solutions while solving smaller systems of equations. The Petra-M code had been used for 3D full wave simulations in minority heating regime for Alcator C-Mod [106,107] and the HHFW regime for both NSTX-U and LAPD plasmas including realistic antenna geometries [108].

VSIM.
Time-domain approaches to solving Maxwell's equations in cold plasma date back to 1994 [109] but suffered from stability issues: the time step had to be chosen small enough to resolve characteristic frequencies in the plasma, such as even the electron plasma frequency or the electron cyclotron frequency, frequencies much higher than those of interest to ICRF. Smithe introduced a modified FDTD scheme which remains stable under the same conditions in plasma as in vacuum [110]. This is achieved by solving Ampère's law together with the constitutive equation for the plasma currents. The two equations are solved implicitly, but their exact discrete locality can be exploited to keep the size of the resulting system of equations constant, not scaling with the total problem size, thus maintaining the desired explicit nature of FDTD. This algorithm has since been implemented in a code called VSIM. The trivial parallelizability of FDTD (thanks to its explicit nature) has enabled the simulation of truly large antennas and large sections of the tokamak [58,111]. The ability, inherent to time-domain approaches, to handle non-linearities, has enabled the inclusion of nonlinear sheath BCs directly [112,113]. An oft-cited disadvantage of the FDTD approach is its reliance on regular Cartesian grids, which prevents resolving find details of a complicated antenna geometry, a task which is much easier with the flexible tetrahedral meshes of finite elements. Attempts to enable more flexible FDTD grids have been somewhat successful [114] but none have been applied to plasmas. Geometrically complex boundaries can be modeled in FDTD using cut-cell techniques [115]. One can formulate the sheath effects as a 'sub-grid sheath BC', in which the effects of the sheath are represented as local capacitive and resistive circuit elements at points on the material boundary [113].
3.1.5. ANTITER. The semi-analytic codes ANTITER II [42] and IV [116] compute in plane geometry the coupling of an antenna array to an inhomogeneous plasma (or dielectric medium). The codes describe the antenna as a set of antenna boxes recessed in the wall, each box being excited by one or several straps. Vacuum is assumed in the boxes. The modelling is made with exp(ik z z + ik y y) Fourier analysis (z and y are respectively toroidal and poloidal directions) and numerical integration along radial direction x. A perfect Faraday screen is assumed at the aperture of each antenna box such that the polarisation of the plasma excitation can be chosen. The antenna array can be tilted with respect to the plasma to take into account the poloidal magnetic field.
ANTITER IV integrates the system of four first order differential equations linking the y and z wave field components. The set of calculated equations describes the fast and SWs with their coupling when their confluence occurs in the density profile. In the limit of vanishing electron mass the set of equations is reduced to the two differential equations that are used in ANTITER II. The reduced equations describe only the FW but they are singular at the Alfven resonance [62] and corresponding to the density region where the confluence with the SW occurs. ANTITER II and IV allow the computation of the strap array impedance matrix, the coupled power to the plasma core and the plasma edge for any current strap distribution. It gives the corresponding k z and k y power spectra and by Fourier inversion the field distribution in the x, y, z coordinates. The two version of codes give about the same results for the coupling, the spectra and even the edge power loss with field-aligned Faraday screen [116]. The codes are relatively fast and can be run on laptop. However, simplification of the 3D antenna geometry to a 2D plane geometry would make the code unable to describe wave excitation as precise as other codes.

Sheath codes
3.2.1. SSWICH. The SSWICH code takes the electric fields as input at the antenna aperture, as calculated by a 3D code without sheaths, such as RAPLICASOL or TOPICA. Then, it re-solves Maxwell's equations in frequency domain, either in full [23] or reduced to only the slow wave component [22], on 2D radial-parallel antenna slices, with the fields from the 3D code at the aperture imposed as BC. The other BCs on the 2D slices are, at first, the asymptotic sheath BC (tangential E fields conservative, D n = 0), and then, as part of an iterative process to determine the sheath thickness δ sh self-consistently, the full sheath BC (tangential E is gradient of V RF ). This is coupled with another part of SSWICH, the DC part, which solves for the DC currents (and the DC potential), excited by the rectified sheath currents, given some anisotropic plasma conductivity sigma ∇ · J DC = −∇ · ↔ σV DC = 0. Work is underway to generalize the SSWICH approach to 3D. The 2D case is fundamentally simpler, since the condition that the tangential electric fields are conservative is a meaningful constraint in 3D but not in 2D. Toy examples have been constructed which show that there are no unsurmountable obstacles preventing a generalization to 3D. It is worth mentioning that SSWICH uses the same multiphysics FEM solver COMSOL as RAPLICASOL, thus they inherently have the advantages of easy meshing and multiphysics coupling, as well as drawbacks of 'black-box' aspects since COMSOL is a commercial software. [29] is a FEM code that is used to explore the interaction of RF plasma waves with material surfaces using a nonlinear sheath model. Such interactions are sensitive to the plasma and RF wave parameters, and the angle that the background magnetic field makes with the material surface. For this reason, rfSOL has been designed with finite element method (FEM) capability to allow flexible surface boundary shapes. ICRF propagation and interaction with the boundary plasma is simulated with a cold fluid model driven by a model antenna with specified antenna current profile. The RF waves couples to the generalized nonlinear sheath boundary including both displacement and particle currents. The rfSOL code is parallelized and uses on the order of 120 cores. It assumes no DC current flow, i.e. the DC sheath BC (equation (13)) is reduced to J n = 0.

Debye scale codes.
In addition, there are a number of codes that model the RF sheath on the Debye scale. These codes are used to obtain the effective RF impedance of the sheath, i.e. its capacitive and resistive properties, which determine the sheath BC for global, device-scale RF modeling. The nonlinear fluid code NoFlu [25], and the particle-in-cell codes hPIC [117] and Vorpal [118] fall into this category. NoFlu solves a 1D nonlinear fluid model for the ions together with Maxwell-Boltzmann electrons to compute the rectified potential and sheath impedance as a function of normalized wave frequency, sheath magnetization, magnetic field contact angle, RF voltage and DC current flow through the sheath. The PIC codes perform the same calculations with kinetic ions (hPIC) and kinetic ions and electrons (Vorpal). Further details are in reference [25]. A body of other work has contributed to the development of micro-scale models and their relationship to global modeling including references [119,120]. [34] is developed for better understanding and modeling of 3D physics in the SOL. It couples the 3D edge plasma fluid code EMC3 [34] and the neutral particle transport code EIRENE [121] self-consistently. EMC3 solves a set of time-independent Braginskii-like equations for mass, parallel momentum, electron and ion energy. EIRENE is a Monte Carlo solver of Boltzmann equation for neutrals. It uses the background plasma parameters calculated by EMC3 as inputs and provides the particle, momentum and energy source terms to EMC3. The parallel transport is treated purely neoclassically in EMC3 while the perpendicular transport coefficients D ⊥ and χ ⊥ are specified as free parameters. These free parameters can be derived through comparisons of simulated profiles with experimental data for current devices and calculated by transport code for future devices. EMC3-EIRENE is particularly robust in calculating 3D inhomogeneous plasma parameters by taking into account arbitrary 3D magnetic field and 3D realistic wall geometries.

Edge plasma fluid/MHD code
Thus, it has been widely used for studying 3D SOL physics in stellarators, tokamaks and linear devices. [54,55] is a non-linear extended MHD code for divertor tokamak studies that is routinely used to investigate large scale plasma instabilities. In particular, ELM and disruption related phenomena as well as the control of these instabilities are addressed with the code. It uses a robust fully implicit time stepping scheme and discretizes plasma and SOL up to the divertor targets via 2D continuous Bezier elements along with a toroidal spectral representation. With this representation, the continuity of the simulation variables as well as their first spatial derivatives is ensured across the computational domain.

JOREK. JOREK
The physics model used in the present study is an energy-conserving ansatz-based reduced MHD model without approximations regarding the geometry [122] with extensions for realistic E × B and diamagnetic background flows [123]. It assumes a stationary toroidal magnetic field and a plasma velocity described solely by parallel flow and E × B plus diamagnetic flows. With such reductions, the model describes time dynamical equations for the poloidal magnetic flux, plasma density and temperature (assuming T e = T i ), electrostatic potential, and plasma velocity in the parallel direction [55]. The implicit time stepping scheme allows for the parallel-to-perpendicular heat diffusion anisotropy to be captured realistically. The JOREK simulation data used in the present study (see section 4) is based on first of a kind realistic type-I ELM cycle simulations [124]. These recover the explosive nature of the instability based on an experimental AUG H-mode scenario.
3.3.3. BOUT++. BOUT++ is a new highly adaptable, object-oriented C++ code for performing parallel plasma fluid simulations with an arbitrary number of equations in 3D curvilinear coordinates using finite-difference methods [56,125]. The goal of the code is to develop a user-friendly, state-of-art, nonlinear fluid turbulence capability for the analysis of SOL turbulence in a general geometry [125]. By separating the complicated (and error-prone) details such as differential geometry and file input/output from the userspecified equations to be solved, the equations being solved are made clear, and can be easily changed with only minimal knowledge of the inner workings of the code. In addition, the efforts of programming can be minimized by calling existing solvers to solve equations. Thus, it allows the user to concentrate more on the physics and develop new physical models. So far, BOUT++ has been successful in modeling turbulence, ELMs, RMPs, peeling-ballooning modes, SOL plasma transport, shear flows, divertor heat flux width in various devices.

ICRF power coupling
ICRF power coupling depends strongly on the width of the evanescent layer and thus the local density in front of the antenna. One of the most effective methods to increase ICRF power coupling is to puff the fueling gas locally close to the antenna. Experimentally, it was found that comparing to the divertor gas puffing, midplane gas puffing close to the antenna increases the coupling resistance by ∼100% (i.e. a factor two) [68]. To reproduce the experiments and understand the corresponding mechanisms, the EMC3-EIRENE code is used to calculate the 3D SOL density, and the antenna code such as FELICE [41] or RAPLICASOL [36] is then used to calculate the coupling resistance with the calculated density.
An example of the calculated density with EMC3-EIRENE for AUG H-mode plasma is shown in figure 2. It is shown that compared to divertor gas puffing (the reference case), top gas puffing increases the density almost toroidal uniformly, but to a small extent. In contrast, midplane gas puffing increases the density very significantly at regions close to the gas valve. The difference in density increase can be explained by the magnetic field line connections. Field lines starting from the top gas cloud toward midplane spread quite widely while the field lines penetrating the midplane gas cloud remain quite concentrated. Consequently, the top gas cloud is connected to all toroidal distributed antennas while the midplane gas cloud is only connected to antennas closest to it. Further studies show that the spreading of the neutral gas, the ionization, the magnetic field line connections to the high density cloud as well as the local decrease of temperature are responsible for the increase of SOL density. The simulated neutral density and electron density were compared with the measured ones and good qualitative/quantitative agreements have been obtained [43,45].
The coupling resistance can then be calculated with the antenna codes such as FELICE or RAPLICASOL. A comparison of experimental and simulated relative increase of coupling resistance (ΔR c /R c_ref ) for both AUG H-mode and Lmode is shown in figure 3. Here R c_ref is the reference coupling resistance during divertor gas puffing and ΔR c = R c − R c_ref . The simulation results are in good qualitative (with 1D FELICE calculations) or quantitative (with 3D RAPLICA-SOL calculations) agreement with the experiments. It is shown that midplane gas puffing close to the antenna increases the coupling resistance most significantly, by ∼100% in H-mode and by ∼60% in L-mode. This coupling resistance increase decreases exponentially as the geometrical distance between the gas valve and antenna becomes larger, and it is related to the transport properties and machine geometries. Thus, it is important to install the gas valve as close to the antenna as possible in the toroidal direction so that the midplane gas puffing can be most effective in increasing ICRF coupling.
Besides local gas puffing, magnetic perturbation (MP) fields can also influence the SOL density and ICRF coupling. MP fields are commonly used in many tokamaks for Comparisons of experimental and simulated relative increase of coupling resistance (ΔR c /R c_ref ) for AUG H-mode and L-mode. MID3 and MID13 represent midplane gas puffing from ports in section 3 and section 13, respectively. Reproduced from [44]. © 2017 EUROfusion Max Planck Institute for Plasma Physics. Reprinted from [68], with the permission of AIP Publishing.
ELM mitigation and suppression [126][127][128][129] and are envisaged for ITER [130]. Their application in H-mode discharges leads to increased particle transport across the separatrix and a non-axisymmetric plasma response [131,132]. Such response kinks the edge flux surfaces in a 3D fashion, creating a 3D density profile that endows the SOL ICRF coupling region with 3D geometry, thus impacting RF coupling characteristics. When MP fields are further rotated, the 3D density profile rotates locked to the external field and coherently affects RF coupling, as first reported in [133]. A dedicated set of discharges to study the effect of MP fields on ICRF coupling was performed in AUG [48]. Here, the resulting 3D density profile in front of the ICRF antennas during the application of MPs could be measured with antenna-embedded X-mode reflectometry. This study led to the conclusion that the measured ICRF evanescent factors (α in R c ∝ R 0 e −αk d evan ) could not be described by well-known formulas [64,70]. Furthermore, fullwave modeling with RAPLICASOL using as input measured 1D density profiles at the midplane could not reproduce the measured coupling resistance in those experiments.
In order to account for the impact of MP-generated 3D geometry on ICRF coupling, a numerical workflow able to reconstruct 3D MHD equilibrium, 3D density and 3D fullwave simulations was devised [47,134]. The PARVMEC and BMW codes were utilized to compute AUG 3D MHD equilibria with applied MP fields. One of these equilibria was used in EMC3-EIRENE to reconstruct the resulting 3D density profile from the separatrix to the vessel wall, which includes the full ICRF coupling region. Full-wave simulations with the RAPLICASOL code were then performed using as input the 3D density profile, finding a remarkable agreement with the measured coupling resistance in the given MP experiment. A Poincaré plot of the PARVMEC & BMW magnetic field with the ECM3-EIRENE grid, as well as a comparison of the modeled 3D density and coupling resistance against experiment can be seen in figure 4.

Slow wave simulations
ICRF antennas do not only excite the FW mode which is used for core plasma heating; they also inadvertently lead to excitation of the SW, either directly or by the confluence with the FW in the SOL. The ICRF SW cannot propagate to high density plasma (roughly above 10 17 m −3 ) and therefore not useful for core heating. By carrying E || which leads to RF sheaths (see section 4.3), damping on collisions and damping close to the LHR [105,116], the SW is among the sources of RF power losses. In most present-day tokamaks there are usually limited possibilities of the SW propagation due to high SOL density. In larger future devices with presumably lower density SOL plasma, the SW losses could play a more important role. Despite the local difficulties with the SW mode below the LH density, the FW coupling to the plasma core remains quite insensitive to the presence of plasma density inside the antenna box, as long as the R cut-off density for the FW at the main k does not move radially [135].
In many simulations of ICRF antennas with codes such as TOPICA or RAPLICASOL, the LH resonance is avoided by including a vacuum layer around an ICRF antenna, and this excludes the region where the SW is propagative. Thus, dedicated numerical studies are required to describe the SW propagation.
Considering the SW launched by an ICRF antenna, the propagating SW forms a 'resonance cone' (RC), a wave with hyperbolic wave-fronts. The wave can be described by the linear wave equation, i.e. by linear combinations of plane waves [62]. It was assessed in [136] that the RC propagation can lead to strong localized voltages and interact with the RF sheaths. For finite elements studies in [137], COM-SOL simulations showed signs of non-convergence even when resolving the small characteristic length scales dictated by the SW dispersion relation (equation (2)). The convergence was achieved by introducing collisions and avoiding vacuum layers for the RAPLICASOL calculations in [105] where a single RC was absorbed before reaching the LH resonance, and in [137] where numerous RCs with multiple reflections were modeled close to the LH resonance (see e.g. [138]). In these cases, the smallest SW length scales were not resolved. Possible reasons for such convergence might invoke strong SW absorption on collisions and an increase of the length scales close to the LH resonance due to local dissipation. The authors of [105,137] hypothesize that it could be sufficient to resolve the characteristic size defined by the scale of variation of the current density on antenna structures. However, it can be argued that the smallest characteristic scale of the SW still needs to be resolved, in order to fully describe the field distribution inside the RC. In addition, sufficient resolution of the (radial) width of the power Figure 5. Radial E-field distribution in a model with propagating FW (above the blue line-FW cut-off) and SW (between the two green lines-SW cut-off and LHR); horizontal cut in the middle of a simplified AUG two-strap ICRF antenna. Reprinted from [137], with the permission of AIP Publishing. absorption peak is important for correct reflection, transmission, and absorption of the wave at the LH resonance, as is addressed in [138]. The correct description of the SW propagation, dissipation and its role in present experiments and future devices is a subject of active research.
In application to AUG, predictions of the SW behaviour are challenging also due to experimental measurements usually targeting densities above 10 17 m −3 , while density data suitable for SW studies is often unavailable and/or unreliable. This can lead to a large range of possible SW field distributions. Coupling the RAPLICASOL SW simulation with codes that provide detailed plasma parameters is beneficial and often necessary. In the example from [137] in figure 5, a profile of electron density from EMC3-EIRENE was employed, as well as the data for neutral particles and ion densities for collisional rates calculation. It allowed obtaining a qualitative picture of both ICRF modes propagating together, without artificial usage of vacuum around the antenna.

ICRF-rectified sheath
Significant E , generated by ICRF antennas at the plasma periphery, is often considered as the source of sheath rectification and enhanced sheath-rectified DC potential. This sheathrectified DC enhanced potential has three main consequences: firstly, the spatial inhomogeneity of the potential can drive E × B drifts across the magnetic field lines and modify the SOL density. Secondly, ions entering the sheath can be accelerated to obtain energy high enough to cause impurity sputtering or self-sputtering. Thirdly, the RF sheath can affect the absorption and reflection of the RF waves at the boundary. The ion acceleration in the RF sheath can result in additional RF power dissipation.
The ICRF-sheath enhanced potential and heat flux as well as the ICRF-sheath driven DC current flows on Toru Supra were simulated with the SSWICH code [22]. The simulation results are in line with the experiments, indicating that SSWICH is able to obtain heat fluxes evolving in opposite directions by considering a DC conductivity tensor and an asymmetry sheath voltage [22,23]. The DC current transport was found as a necessary mechanism to correctly reproduce the radial width of sheath potential even in the presence of FW [23]. With the 3D Monte-Carlo impurity transport and plasmasurface interaction code ERO [139], the RF-sheath enhanced Be erosion on JET outboard limiters connected magnetically to active ICRF antennas was calculated. The results indicate that by adding potentials in the range of 100-200 V (as expected) to the local sheath, the calculated Be sputtering yield matches the measured local Be I and Be II emission, which has an increase by a factor of 2-3 when RF-sheath is present [140]. The impurity sputtering from the plasma facing components of an ICRF antenna can also be calculated by cooperating a fluid plasma model and a plasma-material interaction code Fractal-Tridyn [31]. In addition, a hybrid particle-in-cell model is developed to analyze the kinetic ion energy-angle distributions impacting the RF antenna and its dependence on different plasma and RF sheath parameters [32]. With further benchmarking with experiments, these codes/models can be used as good predictive tools for modeling the RF-sheath and the consequent enhanced potential, heat flux and impurity sputtering.
The AUG three-strap antenna was designed to minimize sheath-rectification-enhanced sputtering, by suppressing the mechanisms which rectify the sheaths, i.e. by preventing the formation of parallel electric fields near the plasma-facing antenna surfaces. Both experimentally and numerically with linear codes, we find that this suppression works best at some specific antenna phasing and power balance [79,141,142]. At these specific settings, the (linear) parallel RF electric field near the plasma-facing antenna surfaces is minimized; the (linear, within the approximations used in the asymptotic [21] SSWICH code [22]) RF potential on those surfaces is minimized; and finally, the nonlinear rectified DC potential, which accelerates ions towards the surface, is minimized.
In [142], we first used the SSWICH code to model all three parts (parallel RF field, RF potential, DC potential) in detail. The results, that is, the optimum operating parameters, are in excellent agreement with experiment. In fact, good agreement is already reached from minimizing the parallel RF field: not much is gained from knowledge of the DC potential in this case. At the time, there were no diagnostics able to measure the DC sheath electric field on AUG. Attempts have since been made to build one based on [143]. When such measurements become available, direct comparison between SSWICH predictions of the DC potential, and measurements, should be possible. Figure 6. Comparisons of the simulated and measured ion saturation current I sat = qn e C s , in which q is the electron charge, n e is the electron density and C s is the ion sound speed. The measured sheath-rectified potential C s is used to calculate the E × B drifts, which are then used in the simulations. Reproduced from [35]. © IOP Publishing Ltd. All rights reserved.
An additional level of sophistication in the modeling of the interactions of RF waves with sheaths is under development, but not yet fully implemented in predictive 3D full waves codes with realistic geometry. In this approach, introduced briefly in section 2, the properties of the RF sheath are represented as a complex-valued surface impedance BC. The effective sheath impedance and rectified potential is determined from a nonlinear micro-scale (i.e. Debye-scale) fluid model which takes into account the profiles of density, electric field, electron, ion and displacement currents flowing across the sheath [24,27]. The resulting sheath BC improves on the previously employed capacitive sheath model, and has been verified by companion PIC modeling codes [25]. When the RF wave solution in the plasma volume is matched self-consistently to the impedance sheath BC, a process which in general requires nonlinear iteration of the RF wave solution with the BC, a more accurate evaluation of the RF wave fields and DC sheath voltage are obtained.
This procedure has been implemented in the 2D code rfSOL [28,29] in simplified geometries. One of the more interesting phenomena that have been seen is the sheath-plasma resonance. This occurs when the sheath capacitance couples to the inductive component of the plasma response in the volume, i.e. the parallel electron current. The RF and DC sheath voltages are significantly enhanced if there is a sheath-plasma resonance. It remains to be seen whether such resonances will be important in realistic antenna geometry. At present, imple-mentation of the full sheath BC into codes such as Petra-M and V-Sim are under development.
Finally, an additional RF sheath related topic that will require significant new code integration efforts is that of global DC current flow driven by RF sheath physics. Previously, local DC current flows near the antenna were explicitly addressed using the simplified anisotropic DC plasma conductivity in the SSWICH code [22,23]. Rectified DC plasma potentials at the sheath boundary in general drive DC currents through the sheath. These currents influence the sheath impedance and the sheath-rectified potential [25,144] depending on their magnitude. The circuit properties must be determined by completing the global current flow around the vessel walls and through the plasma.

ICRF-induced convection
In this section, we will focus on discussing the ICRF induced DC convection with integrated modeling. The interaction between the ICRF waves and SOL plasma is a complex nonlinear process, since the key physical parameters influence each other: the sheath-rectified potential drives drifts in the SOL and changes the SOL density locally; the local density modification will then change the parallel electric field, which will then influence the sheath-rectified potential. To understand this interaction, two numerical methods have been developed: (a) simulations based on experimental measured sheath-rectified potential; (b) self-consistent simulations with several codes, including EMC3-EIRENE [34] (extended to include prescribed drifts [35] and plasma convection), RAPLI-CASOL [36] (ICRF fields) and SSWICH [22] (sheath-rectified potentials).
The first simulation strategy is to use the measurements of the sheath-rectified potential of the reciprocating retarding field analyzer (RFA) to estimate the E × B drifts, recently implemented in EMC3-EIRENE [35]. These drifts are then considered as prescribed parameters and used in EMC3-EIRENE, which was adapted to include prescribed drifts by including the associated drift terms in the mass, momentum and energy equations [35]. This makes it possible to compare the impact of the E × B drifts on the density, predicted by EMC3-EIRENE, with the measurements of RFA or reflectometry. The simulation results are in qualitative agreement with the experimental ones, as shown in figure 6. As pointed out in references [35,145], the results indicate that convective cells are developed where high potential blobs exist. The E × B drifts near the center of convective cells are usually highest (on the order of 1 km s −1 ), and the drifts near the boundary of convective cells can go from one convective cell to the other. The plasma flows driven to the wall by the convective drifts can enhance the interaction between the plasma and wall, and is a possible reason for the enhanced impurity generation and hot spots production on the antenna.
The second method is to iterate the EMC3-EIRENE, RAPLICASOL and SSWICH codes in a self-consistent way to calculate the RF sheath-rectified potential and RF convection. In this method, the 3D density calculated from EMC3-ERIENE is used both in the RAPLICASOL and SSWICH code, the parallel electric field calculated from RAPLICASOL is used in the SSWICH code, and the rectified sheath-rectified potential calculated by SSWICH is used to calculate the E × B drifts, which are then used in the EMC3-EIRENE code as input to calculate the new density. This process is continued until convergence of results is attained using fixed point iteration. An example of this self-consistent simulation loop is shown in figure 7. Three cases have been investigated: (a) twostrap antenna with optimized power ratio P left /P right = 1:1; (b) three-strap antenna with optimized power ratio P center /P outer = 3:2; (c) three-strap antenna with non-optimized power ratio P center /P outer = 1:9. All these studied cases are with dipole phasing. It is shown that the most significant RF convection happens at the top or bottom of the antenna. The two-strap antenna, although with the optimized power ratio between the two straps, produces the largest sheath-rectified potential (>100 V) and RF convection. The three-strap antenna, when operated with the optimized power ratio, generates the lowest sheath-rectified potential and RF convection. This is because with this power ratio, the image currents induced by the center strap and the outer two straps are cancelled by each other [141]. As a result, the image currents on the antenna box are minimized, leading to a minimized parallel electric field and thus a minimized sheath-rectified potential. Figure 8. Mode conversion at a cylindrical density filament. The y-axis of the colorbar shows the normalized parallel electric field (E z ) (normalized to the unperturbed |E ⊥ |) and the x-axis shows to what extent E z is due to the FW or to the SW. The dashed circle line in the figure center represents the filament (radius = 1 cm) boundary. Reprinted from [92], with the permission of AIP Publishing.

ICRF-edge turbulence interaction
In the previous sections, the plasmas are considered at their steady-state and no turbulence is considered. In reality, the SOL plasmas are quite turbulent due to filaments/blobs caused by microinstabilities in any plasma scenario or due to edge localized modes (ELMs) caused by peeling or ballooning unstable modes [146] in high confinement mode (H-mode). The density perturbations caused by filaments or ELMs in the SOL can significantly change the local dielectric properties and influence the propagation of radio frequency (RF) waves. The scattering of RF waves by filaments in the SOL not only influences the amount of RF power flux reaching the absorption region in the plasma core, but may also cause detrimental effects when the scattered RF waves propagate into the plasma facing components. For instance, ICRF waves can induce farfield sheaths and lead to enhanced impurity generation at the wall [147].
Previously, efforts were devoted to understanding the mechanism of RF wave scattering by density filaments [86,90,95,[148][149][150]. Recently, a 3D COMSOL model with a realistic AUG antenna geometry has been built to study the influence of filaments on RF wave propagation [98]. It is benchmarked with the analytical Mie scattering theoretical model [86]. In this model, the filament is aligned with the magnetic field lines and truncated at the boundary by PML. With this model, it is found that the wave scattering depends on the filament radius, density ratio between filament and background density, distance between filaments and number of filaments. In the presence of one filament, a scattering cone is formed. In the presence of multiple filaments, poloidal distributed and radially elongated stripe structures with enhanced and reduced electric fields are developed. Beside the RF wave scattering effect, a filament modifies the power spectrum at regions inside the filament and near the filament-plasma boundary, leading to a redirection of power flow from the perpendicular direction to the parallel direction (parallel to the magnetic field) due to the modification of wave fields inside and near the filament [98]. In addition, significant mode conversion of FW to SW [92,93] can happen near the filament-background plasma interface at resonance conditions, especially when the magnetic field and parallel wavenumber are low.
An example of wave mode conversion is shown in figure 8. The results are generated with an analytic 'Mie scattering' model of ICRF wave-filament interaction. An advantage of the analytic calculation is that the fast and slow modes can be explicitly distinguished, which allows us to show in this figure how SW are generated at the filament interface. The central case (B), where the SW exists in a 'hole'-filament whose density is below that of the surrounding plasma, has been known since 2010 [90]. More detailed calculations revealed that such resonant mode conversion also occurs under other conditions [92,93], i.e. in case A where all waves are evanescent in the filament and background plasma, and in case C where the FW is only propagating in the filament while the SW is evanescent everywhere. This wave mode conversion is now hypothesized to be responsible for several spurious loss phenomena observed in HHFW, helicon, and LH experiments [151].
Furthermore, to understand the influence of ELMs on RF wave scattering in a realistic geometry, the MHD code JOREK is used to calculate the 3D density during an ELM, and the RALICASOL code is then used to calculate the 3D wave fields with the calculated density [53]. An example of such simulations is shown in figure 9. Two cases have been compared, one with a density before the ELM (t = t1) and the other one with a density during the ELM peak (t = t7). The results indicate that a strong perturbation of the SOL density caused by ELM will result in a rather poloidally inhomogeneous wave fields and power flow. Due to the mutual influence of multiple ELM filaments, radially elongated but poloidally distributed stripe structures with enhanced and reduced wave fields develop. Though originating from the SOL, these stripe structures can extend radially into the core plasma and reach the resonance layer. It is expected that they will cause poloidally inhomogeneous heating and change the radial heat deposition profile. Further work are planned in the near future to investigate this hypothesis.
Edge turbulence affects ICRF waves while ICRF can influence edge turbulence in return. A local DC biasing, either due to ICRF-sheath rectification or applied deliberately using an electrode, can influence the SOL turbulence. Experiments in AUG showed that the shear flow generated by ICRF-sheath rectification in the SOL can poloidal stretch or even split the filaments [84]. Simulation with the fluid turbulence code TOKAM has explored the influence of a localized DC bias of the wall on SOL turbulence [152]. The conclusion is that the density inside the convective cells is depleted, the shear flow tears the filaments apart and the fluctuation amplitude decreases locally, which is consistent with the experimental findings. Besides, the turbulence offers a mechanism of transverse DC current transport in the SOL. However one can hardly describe this transport using an Ohm's law. Further studies with a more accurate and realistic model to address this issue is required.

Discussion
In the previous section, modeling strategies and key results of several physical issues of ICRF-edge interaction on AUG are discussed. In brief, they include: (a) ICRF power coupling. This is studied by EMC3-EIRENE + FELICE/RAPLICASOL simulations. With this method, parameters for scenarios with different gas puffing scenarios are calculated, and they are in qualitatively/quantitatively agreement with the experimental data. By validating the simulations with experiments on various devices, the developed tools can be used to find optimized gas puffing locations to maximize ICRF power coupling. Moreover, the PARVMEC and BMW codes are used to compute the 3D induction field in ideal MHD approximation, which is then used in the EMC3-EIRENE and RAPLICASOL codes to study the influence of MP fields on ICRF power coupling. (b) SW propagation. For example, this is studied by the EMC3-EIRENE + RAPLICASOL combination, in which very fine computation mesh has to be built in the vicinity of the antenna (where n e < 10 17 m −3 ) in RAPLI-CASOL to allow proper excitation and propagation of the SW. (c) ICRF-rectified sheath. For instance, this is studied with the self-consistent code SSWICH. Incorporating with the experiments, it has been used to find the optimized antenna feeding scheme on AUG by calculating the DC potential. A more advanced way to calculate DC potential by employing the complex sheath impedance parameter is under development. (d) ICRF induced convection. This is studied with measurements + EMC3-EIRENE simulations or integrated EMC3-EIRENE + RAPLICASOL + SSWICH simulations. With these approaches, the RF enhanced sheath-rectified potential, convective cells and density modification in front of the antenna can be well predicted. Their influence on ICRF power coupling and plasma-wall interactions can also be characterized. Thus, the developed tools from these approaches, and similar ones under development in the USA RF-SciDAC project, can be used to optimize antenna design and plasma scenario in order to minimize RF sheath and plasma-wall interaction. (e) ICRF-edge turbulence interaction. A 3D COMSOL model using realistic power spectrum and experimental density is used to study the power redirection effects by filaments, and an analytical model based on Mie scattering theory is used to study the wave mode convection phenomenon. For a more realistic and comprehensive study, integrated JOREK + RAPLICASOL simulations is performed. With these approaches, the influence of SOL turbulence and ELMs on ICRF wave fields and heating efficiency can be studied.
Among the above mentioned modeling strategies, only the EMC3-EIRENE + RAPLICASOL + SSWICH simulations are run in an iterative way, in which EMC3-EIRENE and SSWICH themselves are iterative and self-consistent codes. Whether an iterative simulation loop is needed will depend on the physics studied. In many cases, coupling an edge plasma code with an antenna code is sufficient, in which the output of the former code is used as input of the latter code. This is because in these cases, the main purpose is to understand how the 3D inhomogeneous SOL density influences the ICRF power coupling or how the 3D time-dependent perturbed density influence the propagation of fast and SWs. The effect of ICRF wave fields on SOL plasma (mainly RF induced SOL convection) is thus neglected. Instead of solving the physical issues independently, it is very important to solve as much physics as possible in one integrated modeling by coupling all necessary codes. Depending on the physics of relevance, the codes available and the difficulties of code coupling, various integrated modeling strategies for ICRF-edge plasma interactions can be developed. An example based on the work discussed in section 4 is shown in figure 10. In this proposed strategy, the physics considered include the antenna physics (full wave excitation and propagation), ICRF sheath and impurity sputtering, ICRF power coupling, ICRF induced drifts and wave-edge turbulence interaction. If considering steady-state plasma, i.e. when plasma turbulence is exempted, the coupling of a 3D edge plasma code, an antenna code and a sheath code is necessary. For instance, these codes could be EMC3-EIRENE + RAPLICASOL + SSWICH. In fact, such selfconsistent simulations were performed previously [38], but the purpose is only to study the ICRF induced convection. It is worth mentioning that EMC3-EIRENE itself is able to simulate impurity transport. However, further development and inclusion of an impurity sputtering model (in the sheath scale) to simulate impurity generation at the wall (especially where the RF sheath-rectified potential is high) in EMC3-EIRENE is necessary. In addition, the RF sheath-rectified potential in front of the wall should be calculated in 3D. By further including the physics of impurity sputtering and transport, the integrated modeling can be made more comprehensive. If considering time-varying plasma parameters, especially when turbulence has to be considered, then the 3D edge plasma fluid code (e.g. EMC3-EIRENE) should be replaced with a 3D turbulence/MHD code (e.g. BOUT++/JOREK) in the simulation loop. Again, impurity sputtering and transport has to be included in the chosen turbulence/MHD code.
The above discussed codes are a non-exhaustive list of codes which have been widely used. They could be replaced by ones which have similar functionalities. In addition, the above mentioned examples are a non-exhaustive list of key physics that could be solved with the existing codes. More physics can indeed be explored with existing or slightly modified codes, for example: (1) reduced impurity production with field-aligned antenna; (2) transport of RF-induced impurity; (3) DC currents due to plasma biasing in presence of mixed metallic and insulating pieces of wall; (4) isotopic effects on RF-induced plasma-wall interaction; (5) density patterns in presence of both RF-induced convection and local gas puffing; (6) assessment of ponderomotive density expulsion in realistic antenna geometry.
Most of the existing codes have been validated with experiments. A few of them, for instance, the sheath, turbulence and impurity transport codes, require further experiments for a more thorough validation. For predictive modeling in future devices, the pre-defined free parameters in some codes, which were determined by experimental profiles or measurements, could be calculated by other code or empirical formula. For instance, the perpendicular transport coefficients in EMC3-EIRENE or other fluid code could be calculated by a transport code. The transverse conductivity in SSWICH could be calculated by an empirical formula which fits the experimental database. On the other hand, these pre-defined free parameters could influence the results quantitatively but not qualitatively. Parameter scans will always be useful for capturing a broader physics picture. To further include ICRF core physics, it would be useful to couple an edge antenna and a core ICRF code in one package, such that the wave excitation, propagation and absorption can all be taken into account. To achieve this goal, the HIS-TORIC (hybrid integration of SOL to TORIC) code package is developed to couple a SOL cold collisional plasma model based on COMSOL with the core heating code TORIC [107,153]. The philosophy is to solve the core and SOL physics separately for all boundary Fourier modes, and then join the two domains together using the continuity BC for the tangential electric and magnetic fields at their interfaces. The inclusion of damping mechanisms in the SOL permits the simultaneous calculations of antenna loading and core heating efficiency. Simulations for Alcator C-Mod show that the computed core RF field pattern with HIS-TORIC is similar to that computed with the regular TORIC code [153]. Later, the Petra-M model [106] (Physics Equation Translator for MFEM) solving the RF wave physics in the SOL was used instead of COMSOL to couple with TORIC [106,107], which allows the simulations to be extended to 3D geometry and run more effectively.
Though the HIS-TORIC or Petra-M-TORIC package calculates the wave propagation and absorption properly, many other physics are neglected. Recently, the Scientific Discovery through Advanced Computing (SciDAC) RF team from USA has made good progress in full integrated modeling [17]. It is aimed to develop simulation capability to predict the selfconsistent interaction of RF power with the core plasma, SOL and wall. As shown in figure 11, this full self-consistent simulation loop mainly couples the SOL and core parts. The SOL part not only has the RF wave solver, but also includes the RF sheath model, the fluid turbulence model, the fluid equilibrium transport model as well as the impurity generation and transport model. The core part contains the RF wave solver and the RF plasma response model. Different models are solved in different dimensions and are allocated with different amounts of Figure 11. Strategy of full integrated modeling including self-consistent interaction of RF power/wave with the core plasma, SOL and wall. This strategy is currently used by the USA RF-SciDAC team. Reproduced with permission from [17]. computational resources. Currently efforts of this project are mainly focused on investigating each subtopic separately by scientists with different expertise. Integrating all related codes or models is still underway.
On the European side, the Integrated Tokamak Modelling Task Force (ITM-TF) is a simulation framework aiming to provide a standardized platform and an integrated modeling suite of validated numerical codes for the simulation and prediction of a complete plasma discharge of an arbitrary tokamak [154]. It uses a data structure called 'Consistent Physical Objects' [155] to couple various physics modules and incorporate both experimental and simulation data. In addition, all machine related data are extracted into standardized machine descriptions so that physics modules can become independent of the specific machine. Sophisticated workflows based on the ITM framework were developed. For instance, the European Transport Simulator [156] workflow uses a sophisticated module to incorporate physical modules of heating, equilibrium, pellets, neutrals, impurities, sawteeth, neoclassical tearing modes and turbulence transport.
The ITER Integrated Modelling & Analysis Suite (IMAS) is a modular set of components enabling collective development and execution of integrated modeling applications [157].
The key component of the IMAS infrastructure is a standardized data model, called the ITER Physics Data Model (PDM). A data access model called application programming interface (API) enables cross-language communication. Physics components can be coupled into an integrated modelling workflow once interfacing to the API data model. Although a large part of the IMAS infrastructure software reuses the ITM developments, significant adaptations and improvements are required to make the original source code consistent with the PDM [155]. Both ITM-TF and IMAS can be considered as frameworks to incorporate multiple RF and plasma codes discussed previously.
Full integrated modeling by coupling multiple codes might be challenging for the following reasons. Firstly, a lot of codes need to be considered in order to account for as much physics as possible. These codes include RF wave codes, RF sheath codes, impurity generation and transport code, kinetic particle or orbit following codes, equilibrium codes, plasma fluid code and turbulence codes. The interfaces between these codes can be very difficult as the physics, grid and time resolution of the codes to be used are quite different. Secondly, the computation resources required may be extremely large, and the convergence of modeling may be rather difficult to reach when a large number of codes are involved.
A more feasible way towards full integrated modeling is to consider each code as a module, in a similar way as the TRANSP code [158]. A 'namelist' can be built to include all input information and have the function of switching on and off the modules. In this way, the computation resources can be best utilized while the computation time can be minimized. For instance, if the purpose is only for studying ICRF power coupling, then only the modules 'EMC3-EIRENE' and 'RAPLICASOL' need to be called. If the influence of ICRF induced drifts on ICRF power coupling need to be further checked, the module 'SSWICH' can be further switched on. In this way, as many codes as possible can be considered, even the same type of code but with different specialty. For instance, 'RAPLICASOL' and 'Petra-M' can be considered as two switches under the same type of code. Nevertheless, multiple interfaces to couple different type of codes need to be deliberately built. Alternatively, multiple codes can be cou-pled into complex simulation workflows utilizing professional tools, with possible switches between codes. For example, both ITM-TF and IMAS employ the KEPLER framework [159].
The other possible approach towards full integrated RF modeling is to develop a very robust code solving a system of equations such that all the above mentioned physics issues are taken into account. These equations mainly include full-wave equations solving wave propagation and absorption, time-dependent fluid and MHD equations solving for turbulent plasma parameters, Fokker-Planck equations solving particle distributions and finally MHD equilibrium equations. The equations for sheath and impurity generation can be considered as BCs, or for more fidelity the RF sheath physics may be also be solved from fundamental equations on a refined mesh. While such treatment of all of these physics problems may make the physics more self-consist and the calculations more efficient, a huge effort may be needed to properly combine all the equations. One of the difficulties is that the timescale considered in different set of equations can have several orders of magnitude difference, and this needs special treatment. The grid also needs special treatment, as the sheath region requires much higher spatial resolution than other part. As for the programming, it can be made much easier with the existing libraries and solvers in MFEM and BOUT++.

Conclusions
Understanding the physics of ICRF-edge plasma interactions is of great importance to the success of ICRF heating in current and future magnetic confinement fusion devices. Because the physics involved is very complex, and very often, different parameters are intrinsically linked, it is only possible to compute them with integrated modeling by coupling several codes, with each code solving a set of equations.
In this paper, the basic physics and modeling codes for ICRF-edge plasma interactions is summarized. The progress in developing integrated modeling approaches to solve several physical issues of ICRF-edge plasma interactions on AUG is then reviewed. They include: EMC3-EIRENE + RAPLICA-SOL to model the excitation and propagation of SW; EMC3-EIRENE + FELICE/RAPLICASOL to model the effect of local gas puffing on ICRF coupling; PARVMEC + BMW + EMC3-EIRENE + RAPLICASOL to model the effect of MP fields on ICRF coupling; RFA measurements + EMC3-EIRENE or EMC3-EIRENE + SSWICH + RAPLICASOL to model the ICRF sheath and convection; JOREK + RAPLI-CASOL to model the influence of ELMs on RF waves. Other analytical methods or standalone numerical codes developed for studying filament induced power redirection and mode conversion are also briefly discussed. Many of the modeling results are compared with the experimental results, and good qualitative and often quantitative agreement is found.
Moreover, possible integrated modeling strategies by including all necessary codes in one package and solving multiple physical issues self-consistently are discussed. For instance, if only considering the physics in the SOL, then by coupling the codes EMC3-EIRENE + RAPLICASOL + SSWICH might be sufficient to solve physics including full wave propagation, ICRF sheath and impurity sputtering, ICRF power coupling and ICRF induced drifts. While impurity transport is considered in EIRENE, a dedicated impurity sputtering model has to be further developed and added in the simulation loop. By replacing EMC3-EIRENE with JOREK, the wave-edge turbulence interaction can be studied. To further include physics in the plasma core, the USA RF-SciDAC team is developing a more comprehensive integrated modeling which has the capability to predict the self-consistent interaction of RF power with the core plasma, SOL and wall. One of the good progresses made within this project is the coupling of the Petra-M (antenna code and SOL wave) and TORIC (core ICRF heating) codes.
The recent progress of modeling ICRF-edge plasma interactions reviewed in this paper is providing better understanding of the key ICRF physics in the SOL. For instance, the integrated modeling to study the effect of local gas puffing on ICRF coupling is playing a very important role in finding optimized local gas puffing locations to maximize ICRF power coupling in current and future fusion machines. With continuing efforts on developing the RF codes and integrated modeling strategies, the modeling tools could be made more robust to solve complex RF physics and design advanced antennas for current and future fusion machines.