Simulation of the ion cyclotron range of frequencies slow wave and the lower hybrid resonance in 3D in RAPLICASOL

Ion cyclotron range of frequencies (ICRF) wave propagation is calculated theoretically for tokamak conditions and for linear magnetized plasma device IShTAR which is dedicated to the RF sheath studies. Only the slow wave (SW) mode of ICRF waves can propagate and be studied in IShTAR. Therefore it is possible to decouple the role of the different ICRF modes in the RF sheath effects. Numerical simulations of the ICRF SW are done in COMSOL in the framework of the existing cold plasma modelling package RAPLICASOL and the SW is for the first time modelled in 3D. To date, RAPLICASOL existed as a 3D wave coupling modelling approach which targets the fast wave (FW). Plasma is implemented as a material with manually assigned physical properties and a perfectly matched layer (PML) is used to absorb the wave energy. Here it is demonstarted how to adjust the RAPLICASOL PML for models with propagating SW. Field structures in the resonance cone shape obtained for the SW differ significantly from the FW and exhibit strong dependence on the density profile in the close proximity of the antenna. The lower-hybrid (LH) resonance is a constant issue in the attempts to model the SW. In this work an approach to obtain correct numerical solutions in the LH resonance presence is demonstrated. Results of this work can be used to improve the complex tokamak ICRF simulations, where so far the SW propagation on the edge has been avoided.


Inroduction
Coupling power to the plasma with ion cyclotron range of frequencies (ICRF) waves is a promising method for heating tokamak plasmas to fusion relevant temperatures. To ensure good wave accessibility to the plasma core, the ICRF antennas must be placed close to the plasma. In experiments on multiple tokamaks, active ICRF antennas modify key plasma edge parameters (plasma density, electron temperature and plasma potential) and enhance destructive plasma-wall interactions. Progress in understanding these interactions is described in details in [1] from the early experiments until the steady conclusion of the RF sheath accelerated ions to be the main reason of impurities release and heat loads.
Direct measurements of the RF-induced parallel electric field would provide quantitative characterization of the RF sheath, but they are hard to perform in a tokamak. One of possible approaches to study parallel electric field of an ICRF antenna directly is to move from a complex tokamak environment to a simple specially designed experiment.
IShTAR (Ion cyclotron Sheath Test Arrangement) [2] is a linear plasma experiment that is close to tokamak conditions for interactions of the ICRF waves and the edge plasma. The ICRF antenna in IShTAR is designed to mimic tokamak antennas, but at the same time its structure is rather simple (no Faraday screen, only one strap). The IShTAR plasma has densities and temperatures similar to the tokamak edge [3]. A stable operational phase (little changes in the plasma parameters such as density and temperature) is significantly longer compared to tokamaks where dramatic changes to the density profile on the edge are introduced on very short time scales. These factors provide favourable conditions for studies of the RF sheath physics on IShTAR.
In order to study the ICRF wave propagation, numerical simulations for conditions typical for IShTAR have been done and are presented here. Theoretical calculations provided below demonstrate that only one of the two ICRF modes, the slow wave (SW), can propagate and be studied in IShTAR. The slow wave propagates at lower densities than the fast wave (FW) and therefore is not used in the present-day big plasma experiments for the purpose of plasma heating or current drive. The existing simulations of the ICRF waves are focused on the FW mode that can deliver power to the hot plasma core. Method of Moments code TOPICA [4] and Finite Element solver RAPLICASOL [5][6][7] evaluate the distribution of the field of the fast wave and introduce artificial modification (vacuum layers) to the density profile in order to avoid the lower hybrid (LH) resonance, which results in the absence of the conditions suitable for the SW propagation. Other simulation tools exist as well, for example an attempt to model the physics of ICRF antenna fields together with plasma-surface interactions using a finite-difference timedomain approach [8]. Lower-Hybrid antennas working in close frequency range are modelled using Finite Element code TORIC [9,10] (which is also used for ICRF modelling) or COMSOL-based LHEAF [11].
However, the SW can play significant role when it propagates in the edge tokamak plasma. The polarization of the FW includes very little parallel electric field component (less than 1 % in the considered example conditions), in contrary to the SW, which can have up to 100% of the parallel component (will be shown below). Therefore the RF sheath on the antenna limiters is usually attributed to the SW. The SW polarization, however, depends strongly on the density and the parallel component is only significant at quite low densities which for tokamaks might be present only in the far limiter shadow region. The relative contribution of the two ICRF modes on the RF sheath formation is unclear; it should vary strongly in the radial direction (because of steep density gradients) and should depend on the ratio of the power carried by each of them. In this work it is attempted to achieve understanding of the SW propagation in order to further study its influence on the RF sheath. In addition, the SW is important because it is able to carry a part of the power emitted by ICRF antennas and this power can be further deposited in unwanted locations of the plasma, instead of its delivery to the plasma core.
The simulations are done in 3D using the simulation tool RAPLICASOL developed in the commercial software COMSOL. The modelling approach presented in this paper consists of 2 main parts: implementation of collisionality and accurate resolution of the specific SW structure. Previous works has encountered difficulties with the SW resonance layer because (a) they did not include sufficient collisionality leading to singular or near-singular fields at the resonance, and (b) because they sought to resolve the 'wavelength' rather than the true spatial scale of the SW resonance cone structure.
The present study is committed to understanding of the SW behaviour and to a global characterization of its field structures and does not include the nonlinear sheath physics.

Wave equation solutions
In homogeneous magnetized cold plasma the wave equation for a plane wave oscillating with angular frequency w as , where k is the wave vector and r is the space vector, is: The dielectric tensor is derived by Stix [12] using the equation of motion of the cold plasma fluid: where e 0 is the permittivity of free space.
The coordinate system in [12], often called the 'Stix frame', is chosen such that B is parallel to z and the two other directions are equivalent for homogenous plasma. It is useful to assume the y-component of the wave vector to be zero, so the perpendicular component is only in x direction and  = +x z n n n . This way we simplify the calculations and this assumption is usually valid for experimental ICRF antennas. The wave equation becomes ( /w = n kc is the index of refraction): The determinant of the left matrix has to be zero for the wave equation to have nontrivial solutions. This requirement gives us the dispersion relation: e  e  e  e e   e   ---------=^^^n   n  n  n  n   n n  n  n  0  4   X   2  2  2  2  2  2   2 2  2  2 which after some transformations looks like: and  e 1 n 2 for both of them. These simplified solutions of the cold plasma dispersion relation are referred to as the fast wave (6) and the slow wave (7).
Within the ideal MHD description there are two types of Alfvén waves: isotropic compressional Alfvén wave and anisotropic shear Alfvén wave. In the wave classification the fast wave belongs to the compressional Alfvén wave type. It means that in the low frequency limit (  w W ci ) the FW is equivalent to compressional oscillations of the magnetic field lines. The slow wave is a shear (or torsional) Alfvén wave. The discussed waves are also called fast and slow magnetosonic waves. The terminology of the wave classification is somehow inconsistent in different literature. The names 'fast wave' and 'slow wave' might be used in a different way from the one defined here. It was, for example, the case in [12] as explained in [13].

Wave propagation: cut-off and resonance
Whenn 2 (ork 2 ) <0, the value ofn is imaginary, so the wave of the form ( · ) k r i exp is evanescent.n 2 has to be positive for the wave to propagate. Transitions between these two states are called resonance and cut-off. A cut-off occurs when = n 0, 2 for the fast wave it corresponds to: And for the slow wave: In experimental devices  n is defined by the ICRF antenna geometry and typically  > n 1. Both the right cut-off (8) and the left cut-off (9) of the FW can be present in plasma. They are usually located in different plasma layers, the right cut-off being closer to the plasma edge. e > 1 is needed for the condition (10) to be achievable, which is only possible when w W > .
ci Such a case might be present in some devices, for example in a tokamak where the magnetic field grows towards the geometrical center of the tokamak, so the cyclotron frequency increases as well and overcomes the ICRF wave frequency at some point. In ASDEX Upgrade tokamak under typical operating conditions the cut-off (10) can be present in the core plasma. Since it happens at quite a large distance from the plasma edge, where the antennas are located, the SW cannot reach this cut-off. The cut-off (11) for the SW lies in the region of small densities; such conditions might be present on the edge of a tokamak as well as in other devices withh low enough densities.
Resonances happen when  ¥ n . 2 For the FW it means  e =n 2 which is equal to (10) and can exist in the tokamak core plasma, as explained above. For the SW a resonance is at and it corresponds to the Lower Hybrid resonance. This resonance happens at higher densities than the cut-off (11).
A common behaviour of the SW in plasma with density gradient and constant magnetic field (figure 1) is that it is evanescent for densities below the cut-off (11), then propagates until the resonance (12) and then it is evanescent again in higher density plasma. The FW often meets the cut-off (9) at densities higher than the SW propagation and transforms at this point from evanescence to propagation. When such a distinction happens, it allows the physics of the propagating SW and the FW to be studied independently in two separate locations. Otherwise both wave modes can propagate simultaneously and then their effects are mixed. Parameters used to plot the figure 1 are typical parameters for ICRF antennas in ASDEX Upgrade tokamak: gas mixture: 98.5% D+1.5% H, frequency: 36.5 MHz, magnetic field: 2 T,  k =8.98 rad m −1 . The magnetic field is considered constant for the simplicity, which leads to the absence of the SW cut-off and the FW resonance defined by (10) as well as the FW left cut-off (9). This simplification does not affect the physics studied in this work.
In figure 1 full solutions of the dispersion relation are also plotted. As we can see in this case, the approximate solutions (6) and (7) match exactly the full roots. This is not always true, but can happen for wide variety of cases.
The y-axis in figure 1 might be not obvious for a reader. The values correspond to the decimal logarithm of | | k 2 multiplied by the sign ofk 2 for each point: This way the part of the plot corresponding to positivek 2 (propagation) lies fully above 0 and the evanescent regions (negativek 2 ) lie below 0. However, (| |) k lg 2 can have negative sign for positivek 2 and then (| |) ( ) <k k lg sign 0. 2 2 * This happens only when < k 1, 2 so for very smallk . 2 We neglect those points and assigne them a value equal to zero (otherwise fork 2 =0.01, for example, the value on y axis would be −2, so a sudden spike would appear). Similarly, all points which have negativek 2 (evanescence) lie only below 0 on this plot. Very few points (sometimes none) are affected by this correction and it does not introduce any inaccuracies. Such an artificial y-axis leads to easily readable visual representation of the data and the mentioned correction is needed to avoid confusing spikes that could be visually mistaken with a resonance, as was the case for example in [14,15].

Electric field polarization
From the eigenvectors (corresponding to eigenvalue 0) of the matrix in (3) it is possible to calculate the relative strength of each electric field component. From the second row of the matrix in the SW approximation: And from the third row of the matrix in the FW approximation: Therefore, the E y component of SW and E z component of FW are negligibly small. In figure 2 an example of eigenvectors values is plotted for the same conditions as in figure 1.
In the range of SW propagation the wave can be considered electrostatic; it satisfies the electrostatic dispersion relation: And since j = -E k i for electrostatic waves, where j is the scalar potential, we get: 1 7 x z Equation (17) defines the direction of the resulting vector of the electric field (and the wave vector) by the ratio of its components. As can be seen from figure 2 and will also be shown below in the plots of the Stix parameters, the electric field of the ICRF FW is purely perpendicular to the background magnetic field with both perpendicular components having similar order of magnitude for the range of densities where the wave propagates. Parallel component E z 2 amounts to less than 0.03% of the total field which makes FW to be a less probable candidate for generation of significant RF sheath along magnetic field lines. It is true at least in the frame of this simplified model, when any realistic non-zero y component of the wave vector produced by antennas in practice is neglected as well as other factors like finite thickness and electrical conductivity of the antenna strap and the presence of various metallic parts inside the antenna box. The simplified analytical model still gives the correct electric field polarization in the region outside of the antenna box on a distance where all local effects become insignificant.
On the contrary, the parallel component of the SW electric field has a non-negligible value in the region of its propagation. At very small densities the parallel direction is even dominant and such field could generate high voltage gradients in the sheath. Then its percentage drops fast with the density growth. In fact, such small densities as 10 13 m −3 are not present in modern tokamaks. The density at the radial position of the antenna limiters is often even higher than the resonance of the SW, so no propagation occurs. Reliable experimental values have not been measured neither for the densities outside of the ICRF antennas in the limiters shadow, nor for the densities in the antenna box. In a tokamak, transient plasma events as ELMs (Edge Localized Modes-shorttime periodic distortion of the plasma boundary), filaments (smaller distortions of the plasma edge), magnetic perturbations, core and edge instabilities constantly change density profiles radially and poloidally during a discharge and happen on different time scales. This is one of the factors that complicate significantly any practical study of RF sheath effects. Densities appropriate for SW propagation can exist at the plasma edge during the discharge; however the exact distribution is unknown.

Shapes of the wavefront
Plasma is one of the materials that can have indefinite permittivity tensor. Such phenomenon is observed in photonic band gap materials and in electromagnetic metamaterials, artificial media created by subwavelength structuring [16,17]. Many unusual properties of these materials are rarely or never observed in nature.
A diagonal matrix is indefinite when one of the diagonal elements has a different sign. The dielectric tensor (2) is not diagonal but is a Hermitian matrix, so the Sylvester's criterion to detemine the definiteness of a matrix can be applied: if all the leading principal minors have positive determinants, a matrix is positive-definite; if the even leading principal minors' determinants are positive and odd ones-negative, a matrix is negative-definite; otherwise a matrix is indefinite.
The elements of the dielectric tensor in the cold plasma approximation depend on the density. The following analysis is performed for the tokamak example conditions that were used in figures 1 and 2. From figure 3 we can see that the determinants of the leading principal minors for (2) have the signs according to the table 1, so e is positive-definite before the SW cut-off, indefinite in the region between the SW cutoff and the SW resonance and negative-definite after (note that in order to characterize the FW correctly a realistic B(r) profile should be taken into account, which is not done here, so the discussion about the FW is omitted). In [16] a 2ndorder equation for the dispertion relation is considered and in this simpler case the definitness of the dielectric tensor points exactly to the form of the equation, which is a hyperbola for indefinite e and an ellipse for definite e. Our case of the 4thorder equation (5) is more complex and the same straightforward relation does not hold true. Generally speaking, in the regions where (2) is indefinite, thek k x z (real) diagram is a combination of elliptic and hyperbolic branches. Luckily, the approximate dispersion relation of the SW (7) is a 2nd-order equation for which it is easy to define the shape from the sign of its discriminant: In the region of the SW propagation > Disc 0, so the equation represents a hyperbola and the wavefront of the propagating SW is hyperbolic. It is still consistent with the fact that the full 4th-order equation solutions are neither purely hyperbolic nor strictly elliptic in general, yet in this particular case it is very close to hyperbolic.
In literature, a hyperbolic wave is often called by a different name-a resonance cone (RC). RCs were first measured in plasma half a century ago [18] and have been studied for many years after (see [19]), also in metamaterials. A specific quality of the RCs is that these waves do not exhibit a periodic structure propagating spatially in the direction of k, on contrary to the more common elliptic waves. The wave energy is concentrated in a thin region of the cone. The field amplitude oscillates inside the cone in space and time, but is negligibly small outside of the cone. The cone thickness is defined by the size of the oscillation source; the size of a cross-section parallel to the magnetic field corresponds to the width of the ICRF antenna strap that excites the wave (will be demonstrated below). The RC size along the wave vector is usually smaller (often much smaller) than the wavelength, so no periodic wave structure is observed.
In plasma, conditions suitable for RCs occur not only in the ICRF range, but also for the Lower-Hybrid range of frequencies. Important modelling of RC interaction with a selfconsistent radio-frequency sheath near conducting walls perpendicular to the magnetic field was done by J.R. Myra and D.A. D'Ippolito [19]. The authors showed that RCs are able to transport strong localized RF voltages to RF sheaths at the wall. When a threshold is exceeded, an order unity fraction of the launched RC voltage couples to the sheath. The threshold is defined analytically in [19] and is a function of the antenna parallel wavenumber, launched RF voltage amplitude, parallel component of the dielectric tensor and plasma parameters. By adjusting all characteristics to stay below the threshold, sheath voltages are expected to be minimized.
ICRF antennas used in hot magnetized plasma devices have different characteristics (frequency, sizes, etc), and so does the plasma in front of them (density, temperature, etc). However, the variations lie in a quite narrow range in respect to the shape of the emitted waves. For most of the magnetic fusion experiments the ICRF waves are elliptic in the region of the FW propagation and hyperbolic (resonance cones) in the region where the SW propagates.

ICRF waves in IShTAR
The IShTAR one-strap ICRF antenna is shown in figure 4. Fed by a coaxial line, a metallic strap couples power inductively to plasma, as in a tokamak, the strap is surrounded by limiters and the geometry is adjusted to the plasma shape to replicate poloidal symmetry along the strap length.
The theoretical description given in the section 2 can be applied to IShTAR in order to investigate the wave behaviour that can be expected in the experiment. A preliminary theoretical study was conducted in [14] before the final design of the ICRF antenna and the achievable range of the plasma conditions were known. Here the calculations are done for the actual present-day IShTAR conditions.
As was mentioned in Chapter 2,  k is defined by the ICRF antenna geometry and typically  > n 1. For a single strap antenna of a width w it is / ~p k w [20], which for w=0.132 m leads to  k theor =23.8 rad m −1 for the IShTAR antenna. However, the spectrum of a real antenna is more complex and it is also affected by the presence of the limiters. The spectrum of  k is found from the COMSOL model of IShTAR ICRF antenna presented in [21]. A 1D profile of the parallel magnetic field component H z is taken in front of the antenna (figure 5) and a Fourier transform is performed ( figure 6). The resulting spectrum has a strong  k =0 component, as expected for a monopole antenna and which is evanescent in our range of densities. This peak has finite width, and one of the points near the peak is  k = 20.93 rad m −1 , the closest to  k , theor is taken for further calculations ofk . 2 The wave propagation for other components of the  k spectrum will be also analyzed. The following conditions are considered for IShTAR: gas -pure Ar or He, frequency-6 MHz, magnetic field-0.1 T.    The full solutions of the wave dispersion relation, as well as the FW and the SW are plotted in figure 7 for argon and figure 8 for helium. The overall behaviour of the two modes is very similar to the tokamak example in section 2, with same cut-offs and resonances present, but at different densities. A yellow box indicates the region of densities accessible on IShTAR, according to the latest experimental data [3]. It can be seen that in IShTAR only the SW propagation is possible, so the RF sheath effects caused only by the SW can be studied, which is a big advantage for the understanding of the physics. The highest possible densities are achieved only in the very center of the plasma. On the edge, near the ICRF antenna, the densities in both Ar and He are below the LH resonance. Further increase of the densities can be possible with the machine upgrade but is not available at the moment. Influence of the change of  k (10.47, 20.93 and 62.79 rad m −1 ) is shown in figure 9 for argon; the tendency is similar for helium. The main difference is the shift of the FW cut-off density value, it is higher for bigger  k . The SW cut-off and resonance remain at the same position. Values of thek are also affected by the change of  k in the region of the SW propagation, important in the present work. The observed tendency gives a valuable conclusion: for FW studies in a device like IShTAR in case of limited range of achievable densities it is beneficial to construct an antenna with lower  k , then it is easier to obtain a region of FW propagation.
If the  k is kept constant, but the magnetic field is varied (0.06, 0.1 and 0.272 T), changes happen not only to the FW cut-off values, but to the SW resonance as well (figure 10). As the magnetic field is lowered, the FW cut-off and the SW resonance move towards each other, until a threshold is reached (B=0.057 T in this case), after which the full roots coincide and the approximate SW and FW solutions become inapplicable (figure 11). For quite an extensive range of the magnetic field (0.1-0.272 T) the differences are very small for the SW in the IShTAR case. For the smaller fields the SW propagation region of densities is extended, but the maximum    density achievable in the experiments anyway lies below the part affected by the magnetic field change.
Electric field polarization, dielectric tensor elements and shapes of the wavefront are qualitatively the same in IShTAR as the results given in the section 2 for the ASDEX Upgrade tokamak example case: the parallel electric field has significant values only for the slow wave and only at very low densities; the propagating SW has the shape of a resonance cone.

3D modelling of the slow wave
A common problem attributed to the numerical simulation of the ICRF slow wave is the task of finding a solution in the vicinity of the Lower-Hybrid resonance due to the solution singularity at this point [15,[22][23][24][25]. As reported in some cases, very short wavelength (because of high perpendicular wavenumber (figures 1, 7, 8)) is a challenge to resolve in a numerical model [15,22]. At least 10 mesh elements per wavelength are needed normally to achieve sufficient convergence of the computations. Calculations on such a fine mesh are impossible with the present-day computers. With a coarser mesh the calculation either does not converge or results in a solution with artifacts that depend on the mesh size, as was shown for example in [15]. In fact, as explained in section 2.4, the characteristic size of a resonance cone wave can be even smaller than the wavelength and then even finer resolution is needed.
The true solution of a wave equation in plasma with a radially varying density needs to be integrable in order to be able to justify the correctness of a numerical solution, which is not the case when a wave propagates through the region of the LH resonance. In this section it will be shown that the specific shape (resonance cone) of the slow wave needs to be considered when discussing a solution near the LH resonance.
It can be possible to obtain a solution in a domain where the SW propagates and the density layer corresponding to the LH resonant is present and such an example will be shown below. It can be done if the propagation region (usually very close to the antenna) is resolved correctly and a collisional term is introduced to exclude the singularity at the resonance. The numerical noise obtained at the attempt to simulate the SW propagation in [15] looks similar to an example that we show below where the RC is not resolved correctly ( figure 16).
It is crucial in the numerical simulations to distinguish clearly the influence of the numerical noise introduced into the obtained solution. Since the problem with non-resolved wave features and associated numerical noise has been observed in the SW simulations before, in this work a stepby-step approach is adopted for increasing the complexity of the model. It will be discussed at each step what numerical noise is observed and whether a solution can be trusted.
The realistic IShTAR geometry with numerous small details was not used for plasma simulations to avoid very high numerical costs and large computational times. Instead a model with simpler geometry was created which retains the main properties of the IShTAR geometry: homogeneity of the plasma, zero wave vector component along y axis and antenna dimensions (and  k ).

Plasma implementation and testing in COMSOL
In the scope of this work COMSOL models have been developed for wave-plasma interaction studies by relying on the main principles of RAPLICASOL [5][6][7]: • Plasma is simulated as a dielectric medium with dielectric tensor derived in the cold plasma approximation; PML was first suggested in 1994 [27] and since then is commonly used to imitate the outgoing propagating wave without any reflection on the interface. It is built-in in the COMSOL software, however the functionality is rather limited and makes it suitable for usage only with standard dielectrics. In RAPLICASOL the PML was for the first time implemented for high-density plasma as a manually defined material [5,6]. The wave equation is solved continuously from the main simulation domain to the PML region, but the coordinate x in the direction of the wave damping is transformed from real to complex stretched coordinate by the stretching function ( ) In a general 3D case where the direction of damping is undefined, the stretching functions are applied as in (18) to all 3 coordinates. The form of stretching functions in RAPLI-CASOL is chosen to be polynomial: PML -the PML depth, p x -the order of the stretching function, ¢ S x -the real stretch and ¢ ¢ S x -the imaginary stretch. A thorough study of the coefficients of the polynomial stretching function and of their influence on the wave reflection from the PML was done in [6]. It allowed to apply in this work the same polynomial function form and to have a starting point for choosing the coefficients.
It was stated in [6] that the PML in the present formalism is not able to handle both forward and backward waves simultaneously. In general, a wave is considered a forward (backward) wave if the scalar product of the phase and group velocities is positive (negative). The fast wave is by definition a forward wave and the slow wave is a backward wave. However, the SW in [6] was taken as an artificially excited single plane wave, so it was not the same as the RC wave considered here. Some more references on the topic of PML for cold plasma can be found in [28].
The presented simulation work shows that the SW in the shape of a resonance cone is not a purely backward wave and in fact can be absorbed by the same PML as the forward elliptic FW in the toroidal direction z (section 4.4). Furthermore, it was observed in our simulations that the SW and FW absorbtion regions are usually spatially separated in the modeled conditions. It will be demonstrated below that for the SW only the toroidal PMLs are crucial, while the radial PML is intended to absorb the FW. It means that both waves can be propagating and coexist in one simulation domain.
COMSOL software provides an opportunity to plot the direction of the Poynting vector and the distribution of the power absorption, so this way the correctness of the chosen stretching function coefficients is double-checked in the performed simulations.
The simplest possible geometry for the model is a perfect electric conductor (PEC) box with an antenna strap (figure 12). The static magnetic field is directed along the z axis and the antenna configuration is typical for tokamaks and for IShTAR: a metallic strap along the poloidal axis (here it is y axis) grounded on the bottom to the box itself and fed with a coaxial line with a port where the wave is excited. The x axis corresponds to the radial direction in toroidal (or cylindrical) plasmas. In this configuration the plasma has a rectangular shape. On the example in figure 12  can be considered. The wave is excited at the coaxial port at the end of the coaxial line. It results in a current distribution in the strap (and the surrounding structures) relevant to the experimental current distribution. It would not be possible in a 2D geometry, since in 2D a current distribution would have to be known in advance in order to calculate a realistic map of the fields created by this current.
In the described RAPLICASOL models plasma is defined as a custom material with the permittivity tensor (2) and the permeability tensor m = 1, according to the theory in section 2 and to [6]. The PML permittivity and permeability correspond to the material which the PML is adjacent to; it was either vacuum or plasma in the performed simulations. Stretching was applied to the z coordinate in the side (toroidal) PMLs, to the x coordinate in the front (radial) PML and to both x and z coordinates in the corner PML regions.
In all following figures of calculated numerical solutions the field E x is plotted (as the strongest or sometimes the only non-zero component), unless other is specified, and all sizes are in meters. In all simulations the density is constant along y and z, it can only vary in x direction (radial).
The model was first tested with parameters typical for RAPLICASOL simulations of the FW. For the conditions characteristic for the ASDEX Upgrade tokamak (gas-D, frequency-36.5 MHz, magnetic field-2 T), the elliptical fast wave mode structure is reproduced successfully in the model as shown in figure 13. The antenna in this model is in vacuum and it can be seen that the propagation starts only in the plasma domain. The density in the plasma is constant  The side view reveals quite good homogeneity of the field along y direction in the plasma ( figure 13(b)). However, it can be seen that periodic spots of the size of the mesh appear on the vacuum-plasma interface. Their size changed when the simulation was performed on a finer mesh and remained equal to the mesh size. In this simulation the spots do not introduce significant disturbance into the propagating wave fields. During the RF phase sweep the spots move vertically along the vacuum-plasma interface but they cause no considerable propagating noise in the x and z directions.
There are 2 types of PML in this simulation: vacuum PML on the sides near the vacuum layer and plasma PML around the plasma domain. The coefficients of the stretching functions were the same in both PMLs and were chosen according to the guidelines in

SW simulations with antenna in vacuum
A defining characteristic of the propagating slow wave is the shape of the wavefront. From the hyperbolic form of the wave equation it follows that the slow wave propagation is restricted to a specific region-the resonance cone (RC), as was explained in section 2. The wave vector k of a resonance cone wave is directed perpendicularly to the cone front. The electric field polarization in this electrostatic wave can be found very simply from (17): E is parallel to k. In figure 14 the antenna is in vacuum and the cone appears in the plasma domain and stretches out at a nearly constant angle, because the density is varied very slowly (linearly) in x direction in a small range: n=6 × 10 11 -10 12 m −3 growing towards positive x. The density assigned in the front PML is equal to the density at the plasma/PML edge. In the side PMLs the density is the same function of x as in the plasma domain.
For n=6 × 10 11 m −3 and  k =20.93 rad m −1 equation (7) yields:k =12.24 rad m −1 , so k should be directed at an angle: 30. 3 20 1 to the magnetic field (z axis). The modelling result in figure 14 shows excellent agreement with this value. For n= 10 12 m −3 :k =23.3 rad m −1 , so a =  45 which is also observed in the plot. If another component of the  k spectrum is taken,k scales accordingly and the resulting angle remains the same.
In figure 15 the RC direction changes much faster with the x coordinate. In this case the density profile is exponential: n=6 × 10 11 -5.4 × 10 16 m −3 . The growing density leads to a rapid increase in thek value ( figure 7), so the vector k shifts towards the perpendicular direction and the cone propagates more to the sides than forward. At the position where the RC enters the side PML the density is n≈2 × 10 13 m −3 . Weak reflections from the side PMLs outer metallic walls appear in this case, and the reflected attenuated wave tends to be closer and closer to the direction parallel to B as it propagates towards higher density. Some noise also originates from the interface between the vacuum and the plasma layers. The reflections should be eliminated in order to get a correct  field distribution in the region of the higher densities. An attempt to improve the PML by making it deeper resulted in a reduction of the amplitude of the reflected wave, but did not eliminate the reflections completely. Further improvement will be achieved at later steps of the modelling when collisions are introduced into the plasma description. However, it can be already seen that, if the refections are not taken into account, the field outside of the resonance cone is considerably smaller than that inside the cone. In the present model the E x near the plasma-front PML interface (at the density near the LH resonance) is around 2 orders of magnitude less than the field in the cone centre. And all other components of the electric field are negligibly small.
In this model the LH resonance layer is present in plasma. However, contrary to previously reported numerical issues [15], no unresolved structures or significant numerical noise is observed in this model. It appears that when the field at the LH resonance is very small compared to the field amplitude in the RC, the LH resonance has no or very little influence on the calculation of the SW field itself. A formally strict case relevant to this situation can be considered: when the field at the LH resonance position is exactly 0. In this case the solution singularity at the LH resonance does not exist and the difficulties associated with it are not present. The presented numerical result does not fulfill strictly the condition = E 0, LH so it can not be claimed that the numerical solution represents exactly the true solution. However, it is an interesting numerical result worth mentioning here. In order to correctly eliminate the LH resonance singularity, we adapt the collisional plasma approach which will be shown below.
The major part of the SW propagation region corresponds to   k k ( figure 7), which means that the wave vector is nearly perpendicular, i.e. the RC is nearly parallel to the magnetic field. In order to see such a structure it is necessary to make a step forward in the model and introduce the strap into the plasma. If a vacuum layer is kept between the strap and the plasma, it is still possible to get a solution in some conditions, but the numerical noise on the vacuumplasma interface, which is negligible for very-low-density plasma (for example 6 × 10 11 m −3 in figures 14 and 15), becomes significant for higher densities. An example is shown in figure 16, where the plasma density on the interface with the vacuum domain is n=10 15 m −3 (the exponential density profile is n=10 15 -1.36 × 10 16 m −3 ). The RC is visible (figures 16(a), (b)) but covered heavily in the other structures which are generated on the border of the vacuum and plasma. It can be seen clearly on the side view ( figure 16(c)) as well. The homogeneity along y is not preserved. The artificial minima and maxima generated on the vacuum-plasma interface (pointed as 'noise' in figure 16(b)) propagate further in the plasma in the x-y direction creating a picture which can be very misleading because it resembles a periodic wave structure. However, no waves are expected to propagate in that region, especially waves of this direction and shape, so it is obvious that the observed structure is a pure numerical noise. The size of the structures attributed to the numerical noise agrees with the mesh size ( figure 16(d)). The RC size in x direction is of the same order as well and it looks poorly resolved.

SW simulations with antenna in plasma
This step forward in the model is achieved by overcoming some difficulties. Since the fields need to be simulated on a very fine scale in x direction and there can be field and/or density gradients along 10 mm of the antenna thickness, the precise geometry of the antenna strap starts to play an essential role, as well as the local mesh on the strap.
A simple strap with a perpendicular cross-section is replaced now by a strap with rounded corners and the edges where the strap is connected to the coaxial feeder are also nicely rounded. The mesh is modified as well; the mesh element minimum size is reduced to 0.8 mm and the maximum size to 10 mm in x direction. The mesh is made nonuniform in the 3 directions. The least important axis is y, along which we expect to see no or very little variations. In the direction of the magnetic field (z) the spatial variations of the fields are also much less significant than in x direction. So the mesh size is multiplied by different factors for the 3 direction: 1-for x, 0.6-for y, 0.75-for z. The final mesh is the densest at the antenna strap and contains at least 7 elements along x on the strap.
A solution with the high plasma density n=5 × 10 15 m −3 (constant in the domain) is demonstrated in figure 17, with the antenna in the plasma. It can be noticed that for this plasma density the RC angle to B is nearly 0. In figure 17 the fields represent 2 effects: the RC originating from the strap and the relatively strong field near the vacuumplasma interface which is an artificial feature due to a high density jump from vacuum to plasma. The high field at the interface is a numerical issue, but it is rather homogeneous along y ( figure 17(b)), so it is not that big of a problem as in the previously considered example of a strap in vacuum. No problem of the artificial periodical structures in the plasma volume occurs, as was the case in figure 16.

Discussion of the PML coefficients
Since it was observed in the shown models that a positive sign of the imaginary stretch coefficients  S x and  S z results in the non-physical negative energy absorption in the front and side PMLs, negative coefficients were used to get positive power dissipation: The following explanation is proposed. For the side PMLs the resonance cone SW is not seen as a backward wave, but as a forward wave because the property of being forward/backward should not be discussed in general, but in relation to a specific direction. It is not the full group and phase velocity vectors that are crucial but only their projections on the axis of stretching which in the toroidal direction result in a forward wave, so the sign of  S z has to be negative (which is also suitable for the FW propagating partially sideways). In the radial direction the SW is a backward wave. However, for the front PML the absorbed power corresponds not to a propagating SW wave (because it does not propagate forward), but to an evanescent FW which is simultaneously excited and carries part of the energy in the forward direction, that is why the sign of the front PML coefficient  S x should be chosen to suit the FW. It is confirmed in a simulation where the plasma domain is replaced with vacuum and no propagation occurs: the absorbed power distribution in the front PML is identical and it can only correspond to an evanescent wave. A comparison is shown in figure 18. The part of energy carried by the FW to the front PML is quite big for models like in figure 17, because the plasma domain has to be narrow in x direction in order to have a fine mesh in the plasma and a long enough propagating region along z. On such a short distance x the FW does not decay much.
The described choice of the PML coefficients made is possible to perform full wave simulations with positive (physically realistic) power absorptions in all PML domains,  which in the past seemed to be not possible because of the limitations on the PML functionality mentioned in [6]. However, in some cases problems can arise. For example, in figure 14 the RC reaches directly the front PML. The coefficients used in the front PML in this model were ( ¢ = S 1, . It resulted in spatially distinct simultaneous negative and positive power absorption in the front PML ( figure 19), where the negative power absorption originated from the FW and a small portion of the positive absorbtion can be clearly localised at the position of the RC (compare figures 19 and 14 for clarity). The FW PML is suitable for the resonance cone SW absorption only in the toroidal direction and is not applicable in the radial direction. While the qualitative picture of the fields looks correct (as discussed in section 4.2), it is not possible to obtain accurate absolute values in such a model because of the negative power absorption (physically equal to power being generated inside PML). Additionally, the RC damping in the front PML looks very inefficient compared to the side PML absorption (figures 18, 19). Luckily, only at very low densities the RC radial propagation is significant. Most cases relevant to experimental ICRF antennas have higher densities and plasma domains large enough radially for the SW RC to propagate sideways (toroidally).
A small portion of negative power absorption is observed also in the side vacuum PMLs (bottom corners in figure 19) and the origin of it is not completely clear. It might originate from the evanescent field going sideways from the strap, but it is not understood why a positive absorption is not achieved. The values of the negative power in these locations are around one order smaller than the values in the other PML domains, so at the moment this contribution is neglected.

SW simulations with collisional losses
What happens in a real device at the location of the side PMLs?If a resonance cone wave meets a metallic wall perpendicular to B, it will be reflected and it will continue to travel until it gets reflected again and again (figure 20). For an RC at a very small angle to the magnetic field the reflection would happen nearly perpendicularly and the cone fronts would overlap with the reflections with very small shifts in x direction. Such wave concentration in a narrow spatial location cannot be resolved numerically. Our attempts to remove the side PMLs and to obtain a solution with the wave reflection from the side walls in relatively high-density plasma (>10 14 m −3 ) failed.
A high-density plasma model without side PMLs cannot converge unless a mechanism of energy absorption is introduced. In the experimental devices the wave energy can be partially transferred to the sheath present on all conducting walls immersed in the plasma. Modelling of this physics was done in [19]. In the current work, the nonlinear sheath physics is not implemented. Another mechanism of the energy absorption is the collisional damping in the whole plasma volume and it is adopted as the next improvement of the discussed model.
The collisions are implemented to the COMSOL model by following replacements (according to the theory in [29]): Corresponding replacements for the ions were omitted, since the collision frequencies of ions with any other particles are negligibly small in the considered case. For electrons an example value of the electron-neutral collision frequency n en =6 MHz was used in all following simulations. In reality the collision frequency should be a function of radius and is not precisely measured in the experimental machines. Our   estimation is that in the limiter proximity it can vary between MHz and GHz range. In this area the neutral-electron collisions should be dominant over others. The chosen value is smaller than p w = f 2 ≈36 × 10 6 rad s −1 , but of a comparable order, so the energy loss through collisions should be significant.
A model without the side PMLs but with the collisional losses was calculated successfully and the result can be seen in figure 21 for the same density n=5 × 10 15 m −3 as in figure 17. The wave decays on the simulated length and no substantial reflections from the metallic walls occur which could prevent the model from converging to a steady-state solution. The absorption of the SW in the plasma volume is clearly seen ( figure 21(b)), as well as the evanescent FW absorption in the front PML and can be compared to the figure 18.

SW simulations with antenna limiters
Since the wave propagates nearly parallel to the magnetic field in the conditions similar to the IShTAR experiments, the wave might never leave the antenna box, since it has limiters on the sides. The limiters can be implemented in the box model in order to obtain the wave behaviour in their presence. The modified geometry is shown in figure 22. The limiters run all the way from the bottom of the box to the top. They are located at the same distance along z from the antenna as in IShTAR and protrude forward as far as the limiters in the experimental device. The corners are rounded to replicate the real limiters and to avoid numerical problems on the sharp edges.
The obtained solution for n=5 × 10 15 m −3 shows that it indeed happens that the wave is concentrated in a very small space in between the antenna limiters ( figures 23(a), (b)). Almost all power is dissipated in that region ( figure 23(c)), the fraction of the FW power reaching the front PML is significantly reduced. For such high density the parallel electric field component E z is very low, one order of magnitude smaller than the main component E x ( figure 24).
For smaller density, for example n=9 × 10 13 m −3 the wave already propagates at an angle and the RC structure is seen outside the antenna limiters ( figure 25). The parasitic structures beyond the RC are of the same size as a mesh element.
From the results of section 4 it can be noticed that the absolute values of the electric field differ a lot in the provided examples. They are dependent on the power coupled to the plasma. In the simulations the power, and as a result the field values as well, have random absolute magnitudes and it was not attempted to adjust them to the real experimental values. The reason is that the experimental values of the coupled power have not been measured yet in IShTAR and it is planned to be done. In general, the generator power is limited to 1 kW. Once the experimental values are obtained, the absolute values in the numerical solutions can be calibrated.

LH resonance in collisional plasma
The same range of densities as in figure 15 with an exponential profile is considered in order to demonstrate a successful SW simulation in the presence of the LH resonance ( figure 26(a)). With collisions included the SW field outside of the cone is negligebly small and at the resonance location (close to the front PML, indicated with the green line) no singularity is present. The RC wave is damped in the plasma volume via collisions and in the side PMLs. The front PML is narrow, since very small evanescent field reaches it (the FW is damped strongly in the thick plasma volume). Note, however, that the PMLs on the sides had to be made wider in z direction (0.4 m against 0.2 m in figure 15) in order to eliminate more efficiently the reflections of the RC from the outer metallic walls (which were present in figure 15).
In figure 26(b) the parallel electric field component is plotted. E z is quite strong (same order of magnitude as E x ) for small densities~10 13 -10 14 where the RC propagates in this model.
To give an idea about the computational resources used, this model has a mesh of 7.2 × 10 5 elements, it was calculated with a direct solver PARDISO on a machine with 24 cores (Each node is equipped with 2 Intel Xeon E5-2680 v3 (2.5 GHz) CPUs) during 45 min, required 280 GB of memory and had 4.7 × 10 6 degrees of freedom.

Conclusions
A complex ICRF waves study for IShTAR is performed. The wave equation solutions are applied to the IShTAR conditions and it is seen that only the SW mode can propagate and be studied in IShTAR, so the present work is fully focused on this mode. The existing modelling package RAPLICASOL for the FW propagation in cold plasma has been extended to include the SW. In IShTAR, same as in the given ASDEX Upgrade tokamak example case, the parallel electric field has significant values only for the slow wave and mostly at very low densities; the propagating SW has the shape of a resonance cone. This proves that the conditions in IShTAR are relevant for studies of the physics of the RF sheath caused by the SW in tokamaks.
In experiments in Ar, densities up to 5 × 10 16 m −3 at the radial location of the antenna limiters are achievable in the present state of the device [3]. For higher densities (>10 14 m −3 ) the wave structure in the experimental device is expected to resemble that of figure 23 (not the absolute values), when the wave does not go outside of the antenna box. The parallel electric field is relatively small in the RC and outside of the box it is nearly zero, so it cannot directly affect the sheath on the outer side of the limiters. Lower densities at the antenna edge and inside the antenna box would result in a resonance cone structure that stretches outside of the antenna limiters ( figure 25). The parallel component of the electric field is stronger at lower densities ( figure 26(b)). Plasma with lower densities can be a more favorable scenario for studies of the RF sheath formation on the outer side of the antenna limiters.
Collisional damping plays a significant role in IShTAR where the neutral pressure is quite high, so the collisional absorption modelled in COMSOL resulted in strong wave amplitude attenuation in the plasma volume and the same can be expected in experiments.
The slow wave propagation was for the first time simulated in 3D. From the obtained RC structures it is obvious that   the simulations of the SW are only possible in 2D or 3D. A 1D model would not be able to represent the RC and to reproduce the wave propagation, in contrast with elliptic waves for which even in 1D a propagating wave can be simulated. The 3D modelling presented here assumed good homogeneity along one coordinate (poloidal). The possibility to control the obtained solutions and to monitor the resulting violations of the assumed homogeneity gave a certainty in the identification of the numerical noise. It is very important that the distinction of the numerical noise in the SW simulations was achieved, since it has always been an issue in previously attempted SW simulations [15,22]. The 3D simulations allow the strap current and the image currents on the limiters to be calculated self-consistently with the fields, which is a crucial difference from a 2D simulation, where the currents have to be known and imposed as initial conditions. A 3D model uses a coaxial port excitation and reproduces the currents and the RCs excited by them realistically in all 3 directions. Therefore, the presented 3D modeling, while implying good poloidal homogeneity, represents important 3D features of experimental antennas, inacessible in 2D simulations.
Critical numerical issues reported in other works in the vicinity of the LH resonance were not observed and some possible numerical noise was demonstrated and explained. A case when the SW propagation is correctly resolved in the antenna vicinity while the LH resonance is present in the simulation domain was shown in figure 26. No significant numerical noise was produced. The characteristic size that needs to be resolved and that defines the wave field distribution is not the wavelength but the antenna radial and toroidal dimensions. By resolving the region of the SW propagation, placing the antenna inside the plasma (instead of additional vacuum layer) and introducing collisions, it became possible to obtain reliable numerical solutions without significant noise. It was also highlighted that the special RC shape of the propagating SW allows a great simplification of the LH resonance handling, since the RC propagates mostly toroidally and does not reach directly the radial location of the LH resonance, which leads to extremely low fields at that location. The crucial precondition for the successful LH resonance treatment is the correctly resolved propagating SW structure.
It was discovered that the same PML can be suitable for absorbing both the SW and the FW and it is implemented as the PML in the toroidal direction. This surprising result contradicts the examples discussed in [6]. It happens because the resonance cone SW, which was not considered in [6] but is excited by realistic ICRF antennas, is a forward wave in the toroidal direction, contrary to the backward SW discussed in [6] which is artificially excited as a single plane wave. The same PML is not applicable to absorb the SW in the radial direction (in that direction it is a backward wave), which is not a big drawback since mostly toroidal propagation of the resonance cone SW is expected in cases relevant to experimental plasma conditions in front of ICRF antennas. Therefore, in the radial direction it is sufficient to place a PML suitable for the FW. The spatial separation of the SW and FW absorption and the correct choice of the PML coefficients described in this paper make it possible to perform full wave modelling with physically correct positive power absorption.
Until present, in the ICRF simulation tools for fusion devices such as RAPLICASOL or TOPICA the edge plasma layer was replaced with a vacuum layer and the slow wave propagation on the edge was avoided. Results of the presented work make it possible to include a propagating slow wave in RAPLICASOL and to improve complex tokamak ICRF simulations.