Self-consistent modelling of radio frequency sheath in
3D with realistic ICRF antennas

L.F. Lu, L. Colas, L. Cao, Guillaume Urbanczyk, B. Lu, Y.Q. Shen, X.J.

Zhang

► To cite this version:

L.F. Lu, L. Colas, L. Cao, Guillaume Urbanczyk, B. Lu, et al. Self-consistent modelling of radio frequency sheath in 3D with realistic ICRF antennas. Nuclear Fusion, 2024, 64 (12), pp.126013. 10.1088/1741-4326/ad7c65. hal-05018593

HAL Id: hal-05018593

https://hal.science/hal-05018593v1

Submitted on 3 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

PAPER • OPEN ACCESS

# Self-consistent modelling of radio frequency sheath in 3D with realistic ICRF antennas

To cite this article: L.F. Lu et al 2024 *Nucl. Fusion* **64** 126013

View the [article online](#) for updates and enhancements.

You may also like

- The geometry of the ICRF-induced wave-SQL interaction. A multi-machine experimental review in view of the ITER operation
  L. Colas, G. Urbanczyk, M. Goniche et al.

- Numerical investigation of a radio frequency sheath with supra-thermal electrons in the presence of the magnetic field
  Jing Ou and Yueheng Huang

- Calculation of the radial electric field with RF sheath boundary conditions in divertor geometry
  B. Gui, T. Y. Xia, X. Q. Xu et al.
# Self-consistent modelling of radio frequency sheath in 3D with realistic ICRF antennas

L.F. Lu<fnref n="1" />*, L. Colas<fnref n="2" />, L. Cao<fnref n="1" />,<fnref n="3" />, G. Urbanczyk<fnref n="4" />, B. Lu<fnref n="1" />, Y.Q. Shen<fnref n="1" /> and X.J. Zhang<fnref n="5" />

E-mail: lulingfeng1988@gmail.com

Received 25 January 2024, revised 5 September 2024

Accepted for publication 18 September 2024

Published 27 September 2024

<figure><img src="image_3.png" /></figure>

CrossMark

## Abstract

Ion cyclotron resonant frequency (ICRF) induced impurity production has raised many concerns since ITER proposed to change the first wall material from beryllium to tungsten. Enhanced DC plasma potential ($V_{DC}$) due to radio frequency (RF) sheath rectification is well known as one of the most important mechanisms behind the RF induced impurities. Our previous work (Lu et al 2018 *Plasma Phys. Control. Fusion* **60** 035003) considered the impact of both the slow wave and the fast wave on the RF sheath rectification in a 2D geometry. It can barely recover the double-hump structure of the $V_{DC}$ poloidal distribution observed in various machines when only the slow wave is modelled using the multi-2D approach which intrinsically assumes the poloidal wavenumber $k_z$ is zero. The fast wave on the other hand is found to be more sensitive to a finite $k_z$ and may need to be tackled in 3D. This work reports our recent progress on the 3D RF sheath modelling. In this new code, the latest RF sheath boundary conditions (Myra 2021 *J. Plasma Phys.* **87** 905870504) and the realistic 3D ICRF antennas are implemented. Compared to the 2D results, the 3D code could well recover the double-hump poloidal distribution of $V_{DC}$ even with the fast wave included, which confirms our speculation on the necessity of treating the fast wave in 3D. While the double-hump pattern is robust in the simulation, the amplitude of $V_{DC}$ is found to be affected by the magnetic tilt angle and the antenna geometry. This emphasizes the importance of adopting a realistic antenna geometry in the RF sheath modelling. The double-hump $V_{DC}$ poloidal structure breaks as the magnetic tilt angle increases. This is explained by the gyrotropic property of the cold plasma dielectric tensor. The spatial proximity effect we identified in the previous 2D simulations is still valid in 3D. Finally, simulation shows the slow wave dominates the RF sheath excitation in the private scrape-off layer (SOL), while the fast wave gradually takes over when moving to the far SOL region. This code could be a new tool to provide numerical support for ITER impurity assessment and ICRF antenna design.

<figure><img src="image_4.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Southwestern Institute of Physics, PO Box 432, Chengdu 610041, China</footnote>

<footnote n="2">² CEA, IRFM, Saint-Paul-Lez-Durance F-13108, France</footnote>

<footnote n="3">³ Dalian University of Technology, Dalian 116024, China</footnote>

<footnote n="4">⁴ IJL, UMR 7198, CNRS-U. de Lorraine, F-54506 Vandoeuvre Cedex, France</footnote>

<footnote n="5">⁵ Institute of Plasma Physics, Chinese Academy of Sciences, Hefei 230031, China</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
Keywords: ICRF, RF sheath, fast wave

(Some figures may appear in colour only in the online journal)

# 1. Introduction

Ion cyclotron resonant frequency (ICRF) is one of the two auxiliary heating methods (the other one is the neutral beam injection) which can heat directly the fuel ions in the present-day tokamaks. Its typical frequency is in the range of 20–70 MHz. Aside from direct ion heating, it has abundant functionalities, e.g., achieving H-mode, driving plasma flow and current, controlling MHD instabilities, mitigating impurities at plasma core, etc. Despite its great importance on magnetic fusion devices, radio frequency (RF) wave interaction with plasma edge often sets limits to the ICRF operation [1, 2]. Recently, concerns of the enhanced impurity excitation have questioned the use of ICRF in ITER after it was proposed to change the main first wall material from beryllium to tungsten. The losses of the ICRF generated fast ions [3–5] and the enhanced DC plasma potential by RF sheath rectification [6, 7] are among the most important physical mechanisms behind the wave-induced impurity generation. One could get hints on the impact of the former mechanism in ITER from the latest JET DTE2 results, which have demonstrated that with a sufficient high plasma density, the fast tritons generated by the 2nd harmonic heating of tritium, which is supposed to be the main ICRF heating scheme for ITER, could be well confined [8]. Several groups from ICRF community around the world are working on the latter mechanism for ITER, i.e., assessing the impurities induced by RF sheath rectification.

The general strategy of modelling the RF sheath rectification includes three processes, i.e., solving the RF field ($E_{RF}$) module, the RF oscillating sheath voltage ($V_{RF}$) module and the DC plasma biasing ($V_{DC}$) module [9, 10]. After knowing $V_{DC}$, one can further evaluate the $E \times B$ convection, the energy of the ions impacting the wall and the DC plasma current. Those quantities are relevant to the density modifications at Scrape-Off Layer region, sputtering yielding of the impurities and the hot spot at the antenna surfaces, respectively [6, 7, 11]. This strategy was first applied in the self-consistent sheaths & waves for ion cyclotron heating (SSWICH) code [9, 10] at CEA. Like many other previous RF sheath codes, which mainly attribute the RF sheath excitation to the parasitically excited slow wave, the first version of SSWICH was in 2D (toroidal-radial) and only considered the slow wave. Sheath was modelled as boundary conditions since its thickness is much smaller than the other scale lengths in the code [9]. Although in 2D, it can reflect the 3D results by scanning independently all the poloidal positions of a ICRF antenna. This method, named as the multi-2D approach, intrinsically assumes that the poloidal wavenumber $k_z$ is zero. After it was developed, the SSWICH-SW (slow wave) version has been tested against measurements on several present machines, including Tore Supra, AUG and JET [9, 12–14]. Later, the fast

wave, which is the main wave excited by the ICRF antenna, was included in the 2D SSWICH-FW (full wave) version [10, 15], it was found that the ‘universal’ double-hump poloidal distribution of $V_{DC}$ observed in various machines, including tokamaks [16–19] and the LAPD linear device [20], cannot be well reproduced by the multi-2D approach when the fast wave is taken into consideration. Moreover, the modelled $V_{DC}$ poloidal structure exhibits strong poloidal modulations, which are not obvious in the experimental observations, e.g. in the field maps measured using a regarding field analyser [9] or an emissive probe [20]. At that time, we conjectured that dropping the finite poloidal wavenumber $k_z$ prevents the multi-2D approach from describing properly the fast wave propagation. Simply extending the sheath boundary conditions (SBCs) used in 2D SSWICH-FW into 3D will lead to the constraints of the normal components of the electric displacement and the magnetic flux density (the so-called DB boundary conditions in electrical engineering [21, 22]) at the sheath boundaries, which is generally difficult to implement, especially for commercial softwares such as COMSOL [23]. To tackle this problem, the Petra-M code avoids directly putting constraints on the normal components by introducing two additional scalars [24]. The other way used by Oak Ridge team is to model the sheath as an additional ‘thin’ layer (although still larger than the real sheath width) on the top of wall material, instead of treating it as boundary conditions [25]. Another merit of a 3D approach with realistic antenna geometry is that it allows addressing simultaneously with the same tool the fast wave coupling and the RF sheath excitation. By contrast, any RF sheath simulation in multi-2D relies on external 3D antenna codes (generally without sheaths) to provide the RF field maps at the antenna aperture.

Clearly different from the above mentioned existing approaches, we used the latest SBCs [26]. It will be shown in the next sections that with the introduction of the plasma current and sheath impedance in the new SBCs, one can avoid the difficulty of constraining both the normal components of the electric displacement and the magnetic flux density, thus paves the way for 3D modelling.

This work shows our latest progress in the development of 3D RF sheath modelling with realistic ICRF antenna structures. The aim of this paper is to take HL-3 parameters as an example to study the parametric dependence of the $V_{DC}$ poloidal structure, the role of the fast wave on the RF sheath excitation, etc. The structure of this paper is the following, section 2 gives a brief introduction of the physics and equations in the 3D code; section 3 demonstrates the capability of this new code to recover the previous simulation results and experimental observations, in particular, the double-hump structure of $V_{DC}$. In section 4, dependences of the $V_{DC}$ structure on antenna geometry and magnetic tilt angle are analysed. The
spatial proximity effect of the fast wave is also discussed. Section 5 revisits the contribution of the fast wave on the RF sheath excitation in 3D. Section 6 is the conclusion.

## 2. General description of the new 3D RF sheath modelling code

Similar to the 2D SSWICH-FW code, the new 3D code solves $E_{RF}$, $V_{RF}$ and $V_{DC}$ modules iteratively based on the COMSOL finite element solver. To ease the numerical convergence, an

asymptotic version is also created. This section describes this new code in details.

### 2.1. RF field module

The RF field module solves the wave equation under the cold plasma approximation,

$$ \nabla \times (\nabla \times \mathbf{E}_{\text{RF}}) - k_0^2 \varepsilon_r(\theta) \mathbf{E}_{\text{RF}} = 0 \quad (1) $$

$$ \varepsilon_r(\theta) = \begin{pmatrix} \varepsilon_{\perp} \sin^2\theta + \varepsilon_{\parallel}\cos^2\theta & i\varepsilon_{\times} \sin\theta & -\varepsilon_{\perp} \sin\theta \cos\theta + \varepsilon_{\parallel} \sin\theta \cos\theta \\ -i\varepsilon_{\times} \sin\theta & \varepsilon_{\perp} & i\varepsilon_{\times} \cos\theta \\ -\varepsilon_{\perp} \sin\theta \cos\theta + \varepsilon_{\parallel} \sin\theta \cos\theta & -i\varepsilon_{\times} \cos\theta & \varepsilon_{\perp} \cos^2\theta + \varepsilon_{\parallel} \sin^2\theta \end{pmatrix} \begin{pmatrix} x_{\text{tor}} \\ y_{\text{rad}} \\ z_{\text{pol}} \end{pmatrix} \quad (2) $$

where $k_0$ is the wavenumber in vacuum. The magnetic field is tilted in the toroidal-poloidal plane as it is in the 2D SSWICH-FW code, see figure 1, where $x$ denotes the toroidal direction, $z$ the poloidal and $y$ the radial direction. $\parallel$ and $\perp$ represent the parallel and the perpendicular directions w.r.t the magnetic field, respectively. The magnitude of the toroidal magnetic field $\mathbf{B}_r$ only depends on the radial position $R$ and it is homogeneous in other directions, i.e. $\mathbf{B}_r = \mathbf{B}_0R_0/R$, where $\mathbf{B}_0$ is the central magnetic field, $R_0$ the major radius. $\theta$ is the magnetic tilt angle in the toroidal-poloidal plane. From figure 1, we have $E_{\parallel} = E_x\cos\theta + E_z\sin\theta$, $E_{\perp} = E_z\cos\theta - E_x\sin\theta$. $\varepsilon_r$ is the cold plasma dielectric tensor under tilted magnetic field shown by equation (2). Note this is expressed in the COMSOL convention, where the RF field oscillates in time as $\exp(j\omega_0 t)$ at the angular frequency $\omega_0$ of the emitted waves. The typical tilt angle for large aspect ratio machines. i.e. WEST and EAST, is 7°, whereas for low aspect ratio machines, i.e. HL-3 (previously known as HL-2M) and AUG, is about 15°.

The geometry of the 3D code has three parts, from bottom to top, the SOL region, the confined plasma region and the perfectly matched layer (PML) region as shown in figure 2. The antenna sits in vacuum and it is surrounded by the SOL plasma. The detailed CAD diagram of the antenna geometry adopted in figure 2 can be found in figure A1. Figure 3(a) gives a zoom in of the SOL region in figure 2, where three pairs of sheath boundaries are considered. The private SOL edges represent the boundaries enclosed by the antenna side limiters. The side limiter edges are the outer peripheries of the antenna side limiters. The far SOL edges correspond to other elements than the antenna that protrude inside the tokamak vacuum vessel, e.g. the main limiters. Similarly, figure 3(b) shows the poloidal sheath boundaries. These sheath boundaries are only enabled when the tilt angle is nonzero. here we have neglected the cross-field sheath [27] due to the difference of the Larmor radius between ions and electrons.

The plasma region in figure 2 consists of the SOL plasma and a part of confined plasma, where a realistic density profile

can be imposed. Since our code only deals with the RF edge physics, no core absorption mechanism is included. Instead, we use the PML technique developed in [28] to artificially damp the wave as if the single pass absorption efficiency of the fast wave is 100%. For the excitation of RF fields, each strap is fed independently by a coaxial line. RF excitation is imposed by setting up voltage or power at the endplates of the coaxial ports. Of course, one can also impose the tangential components of the electric field at the antenna aperture as in the 2D code.

The following two boundary conditions are implemented at the perfect electric conductor (PEC) and sheath boundaries, respectively

$$ E_t = 0 \quad @ \text{PEC} \quad (3) $$

$$ E_t = \nabla_t(z_{\text{sh}}(\theta, V_{\text{DC}}) \cdot J_n(E)) \quad @ \text{sheath} \quad (4) $$

where $t$ represents tangential to the boundary surface. $J_n$ is the plasma current density across the sheath. $z_{\text{sh}}$ is the sheath impedance, which depends on the magnetic tilt angle and $V_{\text{DC}}$,

$$ z_{\text{sh}}(\theta, V_{\text{DC}}) = i \delta_{\text{DS}}(\theta, V_{\text{DC}}) / \omega_0 \varepsilon_0 \quad (5) $$

where $\varepsilon_0$ is the vacuum permittivity. A more elaborate parameterization of the sheath RF impedance is detailed in [29] and can be used in the future works. $\delta_{\text{DS}}$ the width of the Debye sheath, which is given by the Child–Langmuir law under the titled magnetic field,

$$ \delta_{\text{DS}}(\theta, V_{\text{DC}}) = l_{\text{De}}(\theta) \left( \frac{eV_{\text{DC}}}{k_{\text{B}}T_{\text{e}}} + \ln(\cos\theta) \right)^{3/4} \quad (6) $$

In equation (6), $l_{\text{De}}(\theta) = \lambda_{\text{De}}[\cos\theta]^{-1/2}$ is the local Debye length at the edge of the Debye sheath, expressed as a function of the Debye length $\lambda_{\text{De}}$ in the main plasma, $k_{\text{B}}$ the Boltzmann
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Magnetic configuration in the 3D RF sheath modelling code. Adapted from [10]. © IOP Publishing Ltd. All rights reserved.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 2.** 3D geometry of the RF field module. Red dotted line represents the separatrix at the border between the confined plasma region and the SOL plasma region. The bottom plane is the return conductor. Two coaxial ports are left-right aligned. The thin arrow marks the antenna leading edge. Cyan, blue and orange colours indicate the perfectly matched layer, plasma and vacuum region, respectively.</figcaption></figure>

constant, $T_e$ the electron temperature and $e$ the charge unit. Please note that the magnetic tilt angle $\theta$ is defined as the angle between the $B_0$ and the direction $x$, as in figure 1. In the boundary conditions, the same angle or its complementary angle is used as the incidence angle. This is correct because the sheath boundaries are either parallel or perpendicular to $x$. This will not be the case when a curved geometry is introduced.

## 2.2. RF oscillating sheath voltage module

When the RF sheath voltage is large, i.e. $V_{sh} \gg k_B T_e / e$, as it is during ICRF operation, the displacement current is much larger than the conduction current, thus only the former is retained in $J_n$. The RF oscillating sheath voltage is simply obtained from the Ohm' law,

$$V_{RF} = -z_{sh}(\theta, V_{DC}) J_n(E) \quad (7)$$

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 3.** Zoom in of the antenna box and the SOL region of figure 2. Boundaries in blue represent the sheath boundaries. (a) Three pairs of toroidal sheath boundaries, i.e., private SOL edge, side limiter edge and far SOL edge, respectively. Lperp and Lpar are tuned in section 4. (b) Poloidal sheath boundaries, enabled only when the magnetic tilt angle $\theta \neq 0^{\circ}$.</figcaption></figure>

where $J_n = i\omega_0\epsilon_0n \cdot \epsilon_r(\theta) \cdot E$ is the displacement current inside sheath. To save the computation memory, $V_{RF}$ is only solved at the 2D sheath boundaries shown by the blue colours in figure 3.

## 2.3. DC plasma biasing module

The DC plasma potential module solves the conservation of the DC plasma current in the SOL plasma region, as shown in figure 2,

$$\nabla \cdot (-\sigma \nabla V_{DC}) = 0 \quad (8)$$

where the conductivity tensor expressed in COMSOL coordinates (see figure 1) reads,

$$\sigma = \begin{bmatrix} \sigma_{\perp} \sin^2\theta + \sigma_{\parallel} \cos^2\theta & 0 & (\sigma_{\parallel} - \sigma_{\perp}) \sin\theta \cos\theta \\ 0 & \sigma_{\perp} & 0 \\ (\sigma_{\parallel} - \sigma_{\perp}) \sin\theta \cos\theta & 0 & \sigma_{\perp} \cos^2\theta + \sigma_{\parallel} \sin^2\theta \end{bmatrix} \quad (9)$$

where the parallel DC conductivity $\sigma_{\parallel}$ is the Spitzer conductivity and the phenomenological perpendicular DC conductivity $\sigma_{\perp}$ is chosen so that the simulated $V_{DC}$ radial broadening is comparable with the experimental observations.
<table><caption>**Table 1.** 2D SSWICH-FW VS 3D RF sheath modelling.</caption><thead><tr><th>Module</th><th>Physics</th><th>2D</th><th>3D</th></tr></thead><tbody><tr><td rowspan="2">Asymptotic module</td><td>Approximation</td><td>Quasi-insulating</td><td>Quasi-conducting</td></tr><tr><td>Boundary conditions</td><td>DB boundary condition</td><td>Zero tangential electric field</td></tr><tr><td rowspan="4">E<sub>RF</sub> module</td><td>RF excitation</td><td>2D simplified antenna</td><td>3D realistic antenna</td></tr><tr><td>Simulation domain</td><td>2D with k<sub>z</sub> = 0</td><td>3D with finite k<sub>z</sub></td></tr><tr><td>Sheath</td><td>Power lossless</td><td>Joule loss<sup>a</sup></td></tr><tr><td>Sheath boundary condition</td><td>Epol = 0 since k<sub>z</sub> = 0 and E<sub>t</sub> = -∇<sub>t</sub>V<sub>RF</sub>; Constrain of normal D</td><td>Constrain of two tangential E fields by E<sub>t</sub> = -∇<sub>t</sub>V<sub>RF</sub></td></tr><tr><td rowspan="2">V<sub>RF</sub> module</td><td>Simulation domain</td><td>1D curve</td><td>2D plane</td></tr><tr><td>Governing equation</td><td>Line integral of tangential E</td><td>Ohm's law</td></tr><tr><td>V<sub>DC</sub> module</td><td>Simulation domain</td><td>Multi-2D approach to get poloidal distribution</td><td>Direct 3D</td></tr></tbody></table>

<sup>a</sup> Joule loss inside sheath per unit area A due to the presence of RF sheath impedance [26]: $\frac{P_{it}}{A} = \frac{1}{2} \text{Re}(\frac{1}{z_{ah}})|V_{RF}|^2$.

Although the average of $V_{RF}$ along one RF cycle is zero, it can enhance the DC plasma potential via the sheath rectification process. As in the 2D SSWICH-FW code, the sheath rectification is incorporated via the DC SBCs in $V_{DC}$ module,

$$ \left\{ \begin{array}{ll} \mathbf{J} \cdot \mathbf{n} = 0 & \text{Metal} \\ J_{n, \text{DC}} = J_{\text{is}}^{\perp B} \cos \theta \left( 1 - \min \left( \frac{J_{\text{es}}^{\perp B}}{J_{\text{is}}^{\perp B}}, \exp \left( \frac{e(V_f + V_b - V_{\text{DC}})}{K_B T_e} \right) \right) \right) \\ & \perp \text{SBC} \\ J_{n, \text{DC}} = J_{\text{is}}^{\perp B} \sin \theta \left( 1 - \min \left( \frac{J_{\text{es}}^{\perp B}}{J_{\text{is}}^{\perp B}}, \exp \left( \frac{e(V_f + V_b - V_{\text{DC}})}{K_B T_e} \right) \right) \right) \\ & \parallel \text{SBC} \\ V_{\text{DC}} = 0 & \text{Return conductor} \\ V_{\text{DC}} = V_f & \text{Separatrix} \end{array} \right. \quad (10) $$

where $J_{is}^{\perp B}/J_{es}^{\perp B}$ are the ion/electron saturation currents, respectively. The superscript $\perp B$ represents the situation where the wall is perpendicular to $\mathbf{B}_0$. $\perp$ SBC and $\parallel$ SBC refer to the SBC at the boundaries that are perpendicular or parallel to the toroidal direction, i.e., the blue coloured boundaries in figures 3(a) and (b), respectively. $V_f$ is the floating potential and $V_b$ is the extra biasing of $V_{DC}$ due to the presence of RF waves. It is expressed as the modified Bessel function of the first kind $I_0$ [10]

$$ \frac{eV_b}{k_B T} = \ln I_0 \left( \frac{e|V_{RF}|}{k_B T_e} \right). \quad (11) $$

The iterative loop is then closed by updating the sheath width according to equation (6) after obtaining the new $V_{DC}$. These three modules are thus solved iteratively until a convergence is reached.

## 2.4. Asymptotic module

A first guess is needed to ease the convergence of the iterative loop. In contrast to the 2D SSWICH-FW code and also the

Petra-M code [24], where a 'wide sheath' approximation representing the quasi-insulating limit of an RF sheath is used. Either the quasi-insulating limit or the quasi-conducting limit can be used in the asymptotic model of the code to give a first guess for the fully coupled model. The quasi-insulating limit corresponds to a wide sheath regime, where one assumes the sheath width is infinite so as to simplify the RF SBCs $D_n = \frac{e_{sh}V_{RF}}{\delta_{sh}} \to 0$, which leads to the DB boundary conditions in the asymptotic model. On the other hand, the quasi-conducting limit assumes there is no sheath in the asymptotic mode, so that the sheath width is zero. Therefore, from the sheath width point of view, these two asymptotic regimes are opposite. In principle, both of these initial regimes can be used as long as they can get converged and provide a proper first guess. In our case, our code is based on COMSOL finite element method, where implementing the DB boundary conditions turns out to be very tricky. Therefore, we choose the quasi-conducting limit. The asymptotic module solves the asymptotic $E_{RF}$ and $V_{DC}$ in sequence. Since under quasi-conducting limit, all the boundaries are PEC boundaries. $V_{RF}$ is zero so no asymptotic $V_{RF}$ module is needed. Accordingly, the sheath rectification formula in the asymptotic $V_{DC}$ module can be simplified by setting $V_b = 0$ in equation (10). Using the PEC boundary conditions also helps us get rid of the dilemma we had in the 2D SSWICH-FW code where the magnetic tilt angle cannot be exactly zero in order to keep the SBC working correctly [10]. Table 1 lists the difference between the 2D SSWICH-FW code and this new 3D RF sheath modelling code. The physics and the workflow discussed above are summarized in figure 4.

## 2.5. Numerical convergence of the code

The mesh size is 3 mm in the private SOL and 3 cm in free SOL and 2 cm in vacuum. As a comparison, the SW evanescent length $c/\omega_{pe} \sim 2$ cm with density $10^{17} \text{ m}^{-3}$ at the antenna aperture, i.e., at the bottom of the private SOL region. Further refinement to less than 1 cm is conducted in the Faraday screen (FS) bar surroundings and at the sheath boundaries. Imposing further refinement to the mesh does not change the simulation results. The simulations were done in a cluster using 1
<figure><img src="image_9.png" /><figcaption>Figure 4. Physics and workflow of the 3D RF sheath modelling code.</figcaption></figure>

node which contains 80 CPUs and has a maximum RAM of 1024 Gb. E<sub>RF</sub> module requires the majority of memory, approximately 260 Gb. Taking the advantage of COMSOL cluster computing, it only takes about 10 min to solve the E<sub>RF</sub> module. The iterative loop is conducted with the help of COMSOL livelink for Matlab. The convergence is reached when the variation of the variable, i.e. sheath width, between two consecutive iterations is below a prescribed relative tolerance. In COMSOL, let U be the current approximation to the true solution vector and let E be the estimated error in this vector. COMSOL simulation stops the iterations when the pre-assigned relative tolerance, e.g., 0.01 for the fully coupled model, exceeds the relative error computed as the following weighted Euclidean norm [30]:

$$ \text{Relative error} = \sqrt{\frac{1}{M}} \sqrt{\sum_{j=1}^{M} \frac{1}{N_j} \sum_{i=1}^{N_j} \left( \frac{|E_{i,j}|}{W_{i,j}} \right)^2} \quad (12) $$

where M is the number of fields; N<sub>j</sub> is the number of degrees of freedom (DOFs) in field j. The double subscript denotes degree of freedom index (i) and field (j) component. W<sub>i,j</sub> = max(|U<sub>i,j</sub>|, S<sub>j</sub>), where S<sub>j</sub> is the average of |U<sub>i,j</sub>| for all DOFs i for fixed j, times a factor equal to 10<sup>-5</sup> for highly nonlinear problems or 0.1 otherwise [30]. Typically, a convergence can be reached after running a little bit more than 20 iterations with a relative tolerance of 0.01. With the adopted mesh size, running 5 modules once costs about 30 min and convergence can be reached approximately within 8 h after executing about 20 iterations.

### 3. Simulation of RF sheath under HL-3 like plasma parameters: recovery of the V<sub>DC</sub> poloidal structure

The HL-3 tokamak, previously known as HL-2M, was built in 2020 at Southwestern Institute of Physics. The major/minor radius of HL-3 are 1.78 m/0.65 m. The maximum plasma current and magnetic field are 3 MA and 3.0 T, respectively. It has achieved its first 1 MA H mode in August, 2023. In a few years, it will become a key facility to study the burning plasma physics in China. In total about 41 MW auxiliary power including 20 MW NBI, 11 MW ECRH, 4 MW LHCD and 6 MW ICRH are envisaged. Among them, 5 MW NBI, 5 MW ECRH and 2 MW LHCD are presently available.

The ICRF heating schemes and main parameters in HL-3 have been studied in [31]. The RF frequency is around f<sub>0</sub> = 30 MHz for the D<sup>3</sup>He minority heating scheme under B<sub>0</sub> = 3.0 T. Since ICRH is not presently available in HL-3, we will use density profiles obtained from the integrated modelling provided by METIS simulation as shown in [31]. Figure 5 shows the density profile under the baseline scenario with I<sub>p</sub> = 2.5 MA. This scenario is obtained with 2 MW NBI, 3 MW EC and 2 MW ICRH. The last closed flux surface and the antenna leading edge locate at R = 2.4 m and R = 2.43 m, respectively. The magnetic field and plasma density at R = 2.43 m are about 2.2 T and 1.5 × 10<sup>17</sup> m<sup>-3</sup>. A constant extrapolation of this density is used beyond R = 2.43 m. The R cut-off density here is about 3.4 × 10<sup>18</sup> m<sup>-3</sup> for dipole phasing, k<sub>||</sub> = 7 m<sup>-1</sup>, θ = 0° and n<sub>z</sub> = 0, which is also indicated by the black solid line in figure 5. The dielectric tensor component ε<sub>⊥</sub> is about -5 at R = 2.43 m and thus the lower hybrid (LH) resonance (defined as ε<sub>⊥</sub> = 0) [32] is not present in the simulation. In other words, both the fast wave and the
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Edge plasma density profile from METIS simulation for HL-3 discharge. The locations of the antenna leading edge and the R cut-off layer under dipole phasing, $k_{\parallel} = 7 \text{ m}^{-1}$, $\theta = 0^{\circ}$ and $n_z = 0 \text{ m}^{-1}$ are indicated by the vertical lines. A constant density $n_e = 1.5 \times 10^{17} \text{ m}^{-3}$ is imposed beyond R = 2.43 m.</figcaption></figure>

slow wave are evanescent below the R cut-off density and the fast wave becomes propagative above the R cut-off density. Considering a constant density above the LH layer ignores the strong interplay between the fast wave and the LH resonance or even the slow wave cut-off layer, which is important to the local phenomenology. In ASDEX-Upgrade, a density below the LH resonance is attempted to mimic the ICRF specific tungsten sources in ITER [33]. In order to model such a case, one has to take into account the presence of the LH resonance. This will be studied in another dedicated paper.

This section mainly conducts simulations with $\theta = 0^{\circ}$. The RF field is excited by setting the voltage at coaxial port Vport = 1000 V. As a comparison, the RF field gets 10 times larger if one sets the forward power at coaxial port Pport = 1 MW. Both the dipole phasing and the monopole phasing are considered. Figures 6–9 give a general view of the output field maps of the converged 3 modules. The RF field is solved in all 3D domains. Figure 6 is the real part of $E_z$ field at the middle plane, i.e. at the poloidal position z = 0 m, under monopole phasing with $\theta = 0^{\circ}$. $V_{RF}$ is only solved at the 2D sheath boundaries as indicated in figure 3. To ease the interpretation, the CAD diagram of the antenna geometry used in this section as well as the surface current distributions on the strap surfaces are depicted in figure A1. Figure 7 shows the $|V_{RF}|$ at three toroidal sheath boundaries. The field pattern does not change when refining the mesh size. The red localized peaks in the $V_{RF}$ map at left far SOL edge correspond to the locations of the $V_{DC}$ maxima, at the left private SOL edge and the left side limiter edge, sheath plasma waves [34] arise at the sheath boundaries. $V_{DC}$ is solved at 3D in the SOL region. Figure 8 shows a cut of 3D $V_{DC}$ at the middle plane as enclosed by red dashed box in figure 6.

For the moment, the same conductivity values are used as in the 2D SSWICH-FW code, i.e., $\sigma_{\parallel} = 3.5 \times 10^4 \text{ S m}^{-1}$ and $\sigma_{\perp} = 9.3 \times 10^{-2} (2.5 \times 10^{-3}) \text{ S m}^{-1}$ at the private (free) SOL, respectively since the simulated $V_{DC}$ value and its radial broadening is consistent with typical values measured in other machines. For example, the typical $V_{DC}$ values measured by the retarding field analyser in Tore Supra or the emissive probe in LAPD are about several tens of volts to a few hundred volts.

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Real part of the poloidal electric field Re($E_z$) (V m$^{-1}$) at the middle plane, obtained from the converged iterative RF field module. Monopole phasing is used with $\theta = 0^{\circ}$. Excitation is done by imposing 1000 V at each coaxial port. Dimensions of the geometry are also given with a unit of millimetre. Red dashed lines enclose the region where $V_{DC}$ is calculated in figure 8. f = 30 MHz, $B_0 = 3.0$ T, D($^3$He) plasma.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 7.** $|V_{RF}|$ (V) at the left far SOL edge, the left side limiter edge and the left private SOL edge as indicated in figure 3(a), calculated from the converged iterative $V_{RF}$ module. The unit of coordinate is millimetre.</figcaption></figure>

Parallel DC conductivities are so large that the DC potential is homogeneous in the parallel direction. The precise quantitative value of the parallel conductivity does not matter for $V_{DC}$: the plasma behaves as a perfect conductor $|V_{RF}|/B_0$. In the SSWICH-2D code, the perpendicular conductivity is determined in a way that the simulated $V_{DC}$ radial broadening is comparable to the experimental observations in Tore Supra. Whether these conductivity values are suitable for HL-3 tokamak is still an open question. It can only be answered through similar $V_{DC}$ measurement in the future HL-3 experiments as it was done in Tore Supra or other machines.
<figure><img src="image_13.png" /><figcaption>**Figure 8.** Toroidal-radial distribution of $V_{DC}$ (V) in the red dashed box of figure 6, obtained from the converged iterative $V_{DC}$ module.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 9.** 2D (radial-poloidal) maps of $V_{DC}$ (V) under monopole phasing with $\theta = 0^{\circ}$. Converged solution from the iterative model. From left to right: $V_{DC}$ (V) at the left far SOL edge, the left side limiter edge and the left private SOL edge as introduced in figure 3(a). The coordinates are the same as in figure 7.</figcaption></figure>

Figure 8 is very similar to what we obtained previously by the 2D SSWICH-FW code [10], demonstrating the 3D code could well recover the previous results. Figure 9 shows the $V_{DC}$ distribution at the same three toroidal sheath boundaries as in figure 7. One can see the double-hump structure of $V_{DC}$ poloidal distribution which has been documented experimentally in various tokamaks and even in LAPD linear device [20] is now well recovered.

As a comparison, figure 10 shows the $V_{DC}$ distribution obtained by running in sequence the asymptotic model and the iterative model only once. The results in figure 10 are already quite close to the final converged solution in figure 9. Since the converged solution takes relatively a long time, i.e., 8 h with the adopted mesh size, so below we will use the results obtained from running once the asymptotic model and the iterative loop. Moreover, $V_{DC}$ maps from the asymptotic module, where there is no RF sheath rectification present in the calculation, are also shown in figure 11.

An ICRF antenna is designed to excite the fast wave, whose perpendicular electric component is much larger than the parallel one. In monopole phasing with $\theta = 0^{\circ}$, the perpendicular electric field ($E_y/E_z$) is symmetric w.r.t the central axis ($x = 0$) of the antenna, whereas in the dipole phasing, they are anti-symmetric. The amplitude of $V_{DC}$ at sheath boundaries depends on the local electric field. In monopole phasing,

<figure><img src="image_15.png" /><figcaption>**Figure 10.** Same $V_{DC}$ (V) maps as in figure 9, but now solving consecutively the asymptotic model and the iterative loop only once.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 11.** $V_{DC}$ (V) calculated from the asymptotic module, w/o RF sheath, other parameters are the same as figure 10.</figcaption></figure>

the fields excited by the two straps are added up in contrast to partially cancelled out in dipole phasing. Therefore, one can expect that $V_{DC}$ is larger in monopole phasing. Besides, in monopole phasing, due to the lower average $k_{\parallel}$ in the coupled spectrum, the R-cut off layer shifts towards the wall, thereby improving the power coupling. Consequently, for the same RF voltage feeding the antenna, more power is coupled to the plasma under monopole phasing than in the dipole phasing, e.g., the coupling resistance increases from 0.47 Ω to 1.45 Ω when switching from the dipole phasing (figure 12) to the monopole phasing (figure 10) in the simulations. Although the increased coupled power is mostly coupled to the slow wave and may barely contribute to the heating of the plasma at core. We further checked that the $V_{DC}$ amplitudes between dipole phasing and monopole phasing are the same with the same coupled power, as it is found in the experiment that $|V_{DC}|$ is roughly proportional to the square root of the coupled power [35].

Comparing figure 10 with figure 12, one can see that $V_{DC}$ is indeed about 15 V larger in monopole phasing. This is also consistent with the experimental observations [35, 36].
<figure><img src="image_17.png" /><figcaption>**Figure 12.** The same as figure 10, but with dipole phasing imposed at RF feeder.</figcaption></figure>

Note this code does not take into account the realistic power absorption at the plasma core, so both monopole phasing and dipole phasing have the same full power absorption inside the PMLs. Compared to the $V_{DC}$ maps obtained from the multi-2D approach using the 2D SSWICH-FW code, the poloidal modulation often seen at $V_{DC}$ structure has disappeared. This suggests the importance of taking the realistic $k_z$ in the RF sheath modelling especially when the fast wave is taken into account. Note the radial direction in the 3D code directs from plasma to the wall, which is contrary to the 2D code. Thus, the maximum $V_{DC}$ here is in the upper-half plane instead of the lower-half as it is in 2D simulation [9, 10]. Why the maximum $V_{DC}$ and also the maximal plasma-wall interaction have this shape and why they do not locate on the parts closest to the separatrix, but rather near the antenna corners are still badly known. Now with our new 3D code, it is possible to check the parametric dependence of the locations of maximal $V_{DC}$ and get some hints to these questions.

## 4. Dependence of the RF sheath excitation on the antenna geometry and magnetic tilt angle

From our previous work with 2D SSWICH-FW, it is found that both the RF sheath excitation by the fast wave and the slow wave follow the spatial proximity effects [10, 37]: for the slow wave, the RF sheath excitation gets more intense at the sheath boundaries closer to the antenna. For the fast wave, the relevant distance is between the sheath boundary and the fast wave propagation cone. This speculation is backed by the 2D SSWICH-FW simulations: $|V_{RF}|$ at the far SOL sheath boundary decreases with larger toroidal extension, while it increases with larger radial extension as the distance between the sheath boundary and fast wave propagation cone gets smaller. Here we will perform similar exercises to explore the RF sheath excitation at the left far SOL edge with different toroidal and radial dimensions as denoted by Lpar and Lperp in figure 3(a).

Results show that only varying the toroidal and radial extensions generally changes neither the shapes nor the locations of

<figure><img src="image_18.png" /><figcaption>**Figure 13.** $|V_{RF}|$ (V) at the left far SOL edge indicated in figure 3(a) with various toroidal extensions. From left to right, $L_{par} = 1.2$ m, 1.4 m and 1.6 m.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 14.** $|V_{RF}|$ (V) at left far SOL edge indicated in figure 3(a) with various radial extensions. From left to right, $L_{perp} = 40$ mm, 50 mm and 60 mm.</figcaption></figure>

maximal $V_{RF}$ or $V_{DC}$, although it may influence their amplitudes, as illustrated in figures 13 and 14. The maximum $V_{RF}$ appears at the same locations, i.e. $z = \pm 600$ mm, as in the leftmost figure of figure 7 and it does not change with the variations of $L_{par}$ and $L_{perp}$. The amplitude of $V_{RF}$ at the left far SOL edge tends to decrease with the increase of the toroidal extension, whereas it increases with the increase of the radial extension. $V_{RF}$ at other two boundaries has similar behaviours. These behaviours agree with our previous 2D SSWICH-FW simulations. In order to have a better comparison with our previous simulation results, figures 13 and 14 adopt the same Tore Supra like parameters as in [10], i.e., $f = 48$ MHz, $\theta = 0^\circ$, $B_0 = 2.32$ T and $n_e = 10^{18} \text{ m}^{-3}$ at the antenna aperture. The corresponding 2D results are shown in the appendix B. Note quantitative matching between 2D and 3D results is difficult since the two codes adopt different algorithms to evaluate $V_{RF}$ as explained in table 1. Nevertheless, this demonstrates the spatial proximity effect is still valid in 3D and with the full wave excitation, instead of the slow wave only as it is in [37].
We also checked that the double-hump structure still per-
sists when changing the way of RF feeding or the port
alignment (e.g., be it a up-down array or a left-right array).
Figures A2(a) and (b) show the 3D CAD diagrams of another
antenna design for HL-3. The two antenna types have similar
geometric dimensions, as shown by the numbers in figures A1
and A2. For simplicity, like in the previous design presen-
ted by figure A1, all the curvatures and screws are removed
in figure 2 as well. The antenna leading edge locates at
y = -180 mm in the simulation and the radial extension of
the private SOL is 20 mm. In 3D, the coupled power exhibits a
$k_z$ spectrum, unlike in 2D where $k_z = 0$ m⁻¹. Figure A2(c)
also computed $k_z$ and $k_x$ spectra by performing 2D FFT to
the radial coupled power at a x-z cut plane, i.e. (plane with
a constant y) by the following procedures: 1. Extend the pol-
oidal & the toroidal dimensions of the simulation domain; 2.
Define a reference (x, z) plane; 3. Export 2D maps of com-
plex RF field components $E_x, E_z, H_x, H_z$ at the reference plane;
4. Fourier-transform $E_x, E_z, H_x, H_z$ in the two dimensions
(x, z) → (kx, kz); 5. Compute 2D power spectrum for $P_y(k_x,
k_z) = 1/2$ * real($E_xH_z^* - E_zH_x^*$); 6. Integrate spectrum in
$k_x$ or $k_z$ to get a poloidal or toroidal spectrum. The spectrum
result shows the main peak of the toroidal spectrum is about
$k_x = 7$ m⁻¹ and the main peak of the poloidal wavenumber
$k_z$ is indeed at 0 m⁻¹ as it is assumed in 2D. One can see the
spectrum further extends until $|k_z| \sim 20$ m⁻¹. This value is
much smaller than the perpendicular wavenumber of the slow
wave ($k_\perp \sim 374i$ m⁻¹) at the antenna leading edge in the adop-
ted plasma parameter but is comparable with that of the fast
wave ($k_\perp \sim 6i$ m⁻¹) at the antenna leading edge. This suggests
that using a realistic $k_z$ spectrum instead of using $k_z = 0$ m⁻¹
may change the contribution of the fast wave on the RF sheath
excitation.

Simulations with the new antenna geometry suggest that
the modification of the RF feeding structure may also affect
the VDC amplitude at the sheath boundaries. For example,
figure 15 shows the details of the evolution of the connections
between the coaxial line and strap in the RF field module (cor-
responding to the antenna geometry changes from figure A1
to figure A2), the VDC amplitude increases from figure 12
to figure 16. This is probably caused by the reduction of the
strap input impedance, which changes from 0.47 + 43.35i
Ω to 0.16 + 22.28i Ω. On the other hand, the location of
the maximum VDC rather weakly depends on the location of
the RF feeders and the short-circuits. For example, compare
figure 16 with figure 12, we see clearly the poloidal max-
ima are at the same place despite different locations of the
RF feeders and the short-circuits. The reason behind this is
probably due to the fact that these two antenna geometries
have similar current distributions on the strap surfaces, see
figures A1 and A3 for details. Nevertheless, it is crucial to have
a 3D realistic geometry when one wants to have a quantitat-
ive comparison with experiments. It also stresses the import-
ance of taking RF sheath analyses into consideration at the
antenna electrical design stage. RF sheath optimization has yet
to be done for HL-3 ICRF antenna after more realistic plasma
density profiles and optimized antenna geometry (i.e., try to
get larger coupling resistance, lower the maximum voltage at

<figure><img src="image_20.png" /><figcaption>**Figure 15.** Modification of connections between coaxial cable and poloidal current strap from the antenna geometry 1 to the antenna geometry 2. The red panel is the grounded short circuit.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 16.** V<sub>DC</sub> (V) at three toroidal sheath boundaries after the modification of RF feeding scheme as shown in figure 15. The excitation is done by imposing 1000 V at the coaxial ports, dipole phasing, other simulation conditions are the same as figure 12.</figcaption></figure>

feeder and the reflection coefficient, while keeping a modest
$k_{\|} \sim 7-8$ m⁻¹ through tuning the geometric parameters) are
available.

Besides, from equation (10) and also in [38], we know the
RF sheath rectification depends on the incidence angle or the
magnetic tilt angle when the boundary is either perpendicu-
lar or parallel to the toroidal direction as explained in 2.1.
It is interesting to see how the V<sub>DC</sub> poloidal structure shown
in figures 10 and 12 evolves with the variation of the mag-
netic tilt angle. Figure 17 shows with monopole phasing, as
the tilt angle increases, the up-down symmetry of the V<sub>DC</sub>
structure breaks gradually. Here, θ = 7° represents the large
aspect ratio machines, whereas θ = 15° stands for the low
aspect ratio machines. Similar behaviour is found in the dipole
phasing. This is expected due to the factor that the dielectric
tensor equation (2) is gyrotropic w.r.t the toroidal direction
with θ = 0°, i.e., it has rotational invariance along x axis (i.e.,
up-down symmetry), whereas it changes to left-right sym-
metry (rotational invariance along z axis) when θ = 90°.
<figure><img src="image_22.png" /><figcaption>**Figure 17.** Evolution of $V_{\text{DC}}$ (V) poloidal structures at the left side limiter edge of figure 3(a) with variations of the magnetic tilt angle, monopole phasing, other parameters are the same as figure 10.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 18.** Norm of the electric field at the radial position of antenna leading edge, i.e., $y = -180$ mm, with various magnetic tilt angle, (a) $\theta = 0^{\circ}$, (b) $\theta = 7^{\circ}$, (c) $\theta = 15^{\circ}$. Simulation parameters are the same as in figure 17.</figcaption></figure>

Furthermore, in figure 17, the maximal $V_{\text{DC}}$ at the left side limiter edge tends to shift downwards, whereas the maximal $V_{\text{DC}}$ at the right side limiter edge (not shown here) tends to shift upwards. Further examinations of $V_{\text{DC}}$ at the private SOL edge and the far SOL edge show similar behaviours. Moreover, switch from the monopole phasing to the dipole phasing, this behaviour remains the same. Opposite behaviour is observed when reversing the magnetic tilt angle (the orientation of the FS bars also changes), i.e., change from $\theta \to -\theta$. In this case, the maximal $V_{\text{DC}}$ at the left side limiter edge shifts upwards whereas the maximal $V_{\text{DC}}$ at the right side limiter shifts downwards. This shift of maximal $V_{\text{DC}}$ with tilt angle is consistent with the variation of the peak locations of the norm of the electric field at the antenna leading edge ($y = -180$ mm) as marked in figure 2. Figure 18 shows the toroidal-poloidal $|E|$ map at the radial location of the antenna leading edge by solving once the asymptotic module and the iterative RF field module. Numbers indicate the geometric dimensions. One can see with the increase of the magnetic tilt angle, $|E|$ at the left boundary of the side limiter indeed shifts downwards, while $|E|$ at the right boundary of the side limiter shifts upwards.

### 5. Contributions of the perpendicular and parallel electric fields on the RF sheath excitation

RF sheath excitation is often attributed to the RF field component $E_\|$ in the traditional slow-wave-only RF sheath modelling. To have a quick view of the impact of the $E_\|$, we compared $V_{\text{DC}}$ structure with and without FS. Results show that the maximum amplitude of $V_{\text{DC}}$ increases about 20% at the private SOL edges after removing the FS bars, whereas it has little impact on the $V_{\text{DC}}$ at the side limiter edge and the far sol edge shown in figure 3(a). This is a first indication that the $E_\|$ mainly affects the RF sheath excitation in the private SOL region. A further test is done by disentangling the $E_\|$ and $E_\perp$ components at the excitation by imposing only the corresponding electric field component at the antenna aperture, i.e. at the bottom of the private SOL region in figure 3(a). Figures 19 and 20 show the $V_{\text{DC}}$ structure at three sheath boundaries with only $E_\|$ and $E_\perp$ field map imposed at the antenna aperture, respectively. One can see the up-down symmetry of $V_{\text{DC}}$ poloidal structure appeared in figure 10 deteriorates when there is only $E_\|$ or only $E_\perp$ excitation.
<figure><img src="image_24.png" /><figcaption>**Figure 19.** Same as figure 10, but with only $E_\parallel$ field excitation at the antenna aperture.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 20.** Same as figure 10, but with only $E_\perp$ field excitation at the antenna aperture.</figcaption></figure>

Figure 19 is very similar to what we obtained from the 2D SSWICH-SW code. Compared to figure 10, the maximal $V_{DC}$ at the private SOL edge with $E_\parallel$ only excitation in figure 19 just decreases by about 20 V, whereas in figure 20, it decreases by about 40 V when imposing only $E_\perp$. This further confirms that $E_\parallel$ dominates the RF sheath rectification in the private SOL, whereas the $E_\parallel$ has less impact on the RF sheath excitations on the side limiter edge and the far SOL edge: from figure 19 $V_{DC}$ only increases by 15 V compared to the case without RF sheath as in figure 11. On the contrary, $V_{DC}$ is larger under $E_\perp$ only excitation at the far SOL edge and the side limiter edge, suggesting that $E_\perp$ may play a major role on the RF sheath excitation at these two boundaries. This is consistent with our expectation that the slow wave which has a major $E_\parallel$ component dominates the RF sheath excitation at the antenna near field, whereas in the far SOL region, the fast wave becomes the governing wave since the slow wave is quickly vanishing away from the antenna and cannot propagate to the far SOL region.

Another crude way to confirm the slow wave having the dominant role on the RF sheath excitation at the private SOL edge is to artificially modify the dielectric tensor elements $\epsilon_\times$ and $\epsilon_\parallel$. The evanescence length of the slow wave is inversely proportional to $\sqrt{|\epsilon_\parallel|}$, so reducing $\epsilon_\parallel$ will make the slow wave more propagative while it has almost no impact on the propagation of the fast wave due to the distinct dispersion properties between these two waves. However, one must be aware of the mutual coupling between these two waves at the tilted magnetic field configuration [39] or when there is a sharp transition of the geometry [40]. Simulation shows the amplitude of $V_{RF}$ at the private SOL edge indeed increases with the reduction of $\epsilon_\parallel$, while it remains almost invariant with the reduction of $\epsilon_\times$.

In summary, our results in this section suggest the underlying assumption made in the previous slow-wave-only RF sheath simulations that neglecting the contribution of the fast wave on the RF sheath excitation may be valid only when dealing with the private SOL. When it goes to surfaces
hardly accessible to the evanescent slow mode, starting from
the outer edges of the ICRF antenna side limiters, one
should already consider simultaneously the fast wave in the
calculation.

6. Conclusion

This paper reports our recent progress on modelling the RF
sheath in 3D. After applying the latest RF SBCs and the real-
istic 3D ICRF antenna geometry, the V<sub>DC</sub> poloidal modu-
lation observed in the previous 2D SSWICH-FW code has
been eliminated. Three-dimensional (3D) full wave simulation
has well recovered the V<sub>DC</sub> poloidal structure that has been
observed in many machines. Besides, the dependence of this
poloidal structure on the magnetic tilt angle and the antenna
geometry was also explored. Results show that the distribu-
tion of V<sub>DC</sub> is sensitive to the magnetic tilt angle but weakly
depends on the way of connecting the feeding coaxial cable
to the strap or the locations of the RF feeders and the short-
cuts, although the latter would affect the amplitude of V<sub>DC</sub>
through modifying the strap input impedance. This suggests
that modelling RF sheath in 3D with realistic magnetic config-
uration and antenna geometry is not only legitimately import-
ant, but also a must especially when taking the fast wave into
consideration.

One of the interesting features of the RF sheath excitation
is the spatial proximity effect, which has been identified previ-
ously by our 2D simulation. It is found that this effect persists
in our 3D results. Regarding the contribution of the fast wave
on the RF sheath excitation, simulation confirms that in the
private SOL, the slow wave is indeed the dominant wave that
contributes to the amplitude of the RF sheath voltage, which
validates the methodology adopted by most of the previous RF
sheath codes that model only the E//. As for the side limiter
edge or the far SOL region, the fast wave must be taken

into consideration. However the RF sheath oscillating
voltages on these surfaces are generally of lower magnitude
than in the private SOL. In addition, the surfaces connecting
magnetically to the antenna limiters may ‘feel’ the DC plasma
bias induced in the private SOL via the DC current transport.

This code may help the HL-3 ICRF antenna design, in par-
ticular, to minimize the RF sheath by optimizing the antenna
geometry when a more realistic density profile that is repres-
entative to the plasma performance with ICRF on is avail-
able. Also, this 3D code needs more validations by making
simulations that could be compared with other codes [41]
and experiments. For future work, it would be interesting to
make simulations that could be compared to HL-3 ICRF exper-
iments. This is the best validation of the simulation tools.
Moreover, it would be also worth repeating the exercise with
surfaces at shallow angle with respect to *B*, and in the pres-
ence of propagative slow mode. Although this code uses HL-3
parameters in the analyses, it can be easily applied to other
machines. Works have been done using the 2D SSWICH-SW
code by the ICRF community to assess the impurity sputtering
for ITER under the tungsten first wall [42]. Now with this new
code, it is possible to remove some approximations made in
the previous study and strengthen our confidence in the RF-
sheath estimates.

Acknowledgments

This work is supported by the National Magnetic Confinement
Fusion Science program of China under the Contract Number
2022YFE03090002, the National natural Science Foundation
of China under Grant Agreement 12375212, China National
Nuclear Corporation Fundamental Research Program under
Contract Number CNNC-JCYJ-202236 and the Sichuan
Science and Technology Program with Contract Number
2023ZYD0016.
Appendix A. CAD diagrams of the antenna geometries used in this paper

<figure><img src="image_26.png" /><figcaption>**Figure A1.** (a) and (b) Front and lateral views of the antenna geometry 1 used in figure 6. The dimensions are marked in millimetre. The blue boundaries are the strap front surfaces where the surface currents are plotted in (c). Red panels are the grounded short circuits. The private SOL region is also indicated in (b). (c) Surface current distribution under monopole phasing, other plasma parameters are the same as in figure 6.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>**Figure A2.** Schematic view of the ICRF antenna geometry 2 used in figure 16 and the expected power spectra. (a)–(b) Back and lateral views of the two strap antenna. (c) Poloidal and toroidal power spectra, dipole strap phasing, *f* = 30 MHz, *B*₀ = 3.0 T, D(³He) plasma.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Figure A3.** Left: front view of the antenna geometry 2. Blue boundaries are the strap front surfaces where the current distributions are plotted on the right figure. A zoom in of the feeder connection to the strap is shown. To have a better visualization, the antenna septum has been hidden. Right: Surface current distribution under dipole phasing, other parameters are the same as in figure A2.</figcaption></figure>
## Appendix B. Demonstration of the fast wave proximity effect in 2D

<figure><img src="image_29.png" /><figcaption>**Figure A4.** (a) Geometric dimensions of the 2D SSWICH-FW code. The background field pattern is the poloidal electric field with Poynting vector (red arrow). The boundary enclosed by the black ellipse is the target sheath boundary where $V_{RF}$ is picked. (b) Variation of $|V_{RF}|$ with different Lpar. (c) Variation of $|V_{RF}|$ with different Lperp. (b)–(c) Use the same plasma parameters as in figures 13 and 14.</figcaption></figure>

## ORCID iDs

L.F. Lu https://orcid.org/0000-0002-2710-8826

L. Colas https://orcid.org/0000-0002-4573-3326

G. Urbanczyk https://orcid.org/0000-0002-4026-6802

## References

[1] Hillairet J. 2023 *Rev. Mod. Plasma Phys.* **7** 16

[2] Noterdaeme J.M. 2020 *AIP Conf. Proc.* **2254** 020001

[3] Zheng Y., Zhang W., Li G., Yin L. and Hao B. 2023 *Nucl. Fusion* **63** 046016

[4] Mantsinen M.J., Graves J.P. and Cooper W.A. 2021 *Plasma Phys. Control. Fusion* **63** 094002

[5] Kazakov Y. et al 2017 *Nat. Phys.* **13** 973

[6] Colas L. et al (JET Contributors, the ASDEX Upgrade Team, the EAST Team, the WEST Team and ITPA IOS) 2022 *Nucl. Fusion* **62** 016014

[7] Urbanczyk G., Fedorczak N., Gunn J., Colas L., Li J.G. and Wang K. 2021 *Nucl. Mat. Energy* **26** 100925

[8] Mantsinen M.J. et al (JET Contributors) 2023 *Nucl. Fusion* **63** 112015

[9] Jacquot J., Milanesio D., Colas L., Corre Y., Goniche M., Gunn J., Heuraux S. and Kubič M. 2014 *Phys. Plasmas* **21** 061509

[10] Lu L. et al 2018 *Plasma Phys. Control. Fusion* **60** 035003

[11] Beers C.J. et al 2021 *Phys. Plasmas* **28** 103508

[12] Zhang W. et al 2022 *Nucl. Fusion* **62** 075001

[13] Bobkov V. et al 2019 *Nucl. Mat. Energy* **18** 131

[14] Krivska A., Bobkov V., Colas L., Dumortier P., Durodié F., Jacquet P., Klepper C.C., Milanesio D. and Urbanczyk G. 2019 *Nucl. Mat. Energy* **19** 324

[15] Colas L. et al 2017 *EPJ Web Conf. vol* **157** p 01001

[16] Colas L. et al 2007 *Plasma Phys. Control. Fusion* **49** B35

[17] Bobkov V. et al (ASDEX Upgrade Team) 2010 *Nucl. Fusion* **50** 035004

[18] Klepper C.C. et al 2013 *J. Nucl. Mater.* **438** S594

[19] Urbanczyk G. et al 2020 *Nucl. Fusion* **60** 126003

[20] Martin M.J., Gekelman W., Van Compernolle B., Pribyl P. and Carter T. 2017 *Phys. Rev. Lett.* **119** 205002

[21] Lindell I.V. and Sihvola A.H. 2009 *IEEE Trans. Antennas Propag.* **57** 694

[22] Lindell I.V. and Sihvola A.H. 2010 *IEEE Trans. Microw. Theory Tech.* **58** 903

[23] Lu L. 2016 Modelling of plasma-antenna coupling and non-linear radio frequency wave-plasma-wall interactions in the magnetized plasma device under ion cyclotron range of frequencies *PhD Thesis* University of Lorraine

[24] Shiraiwa S., Bertelli N., Tierens W., Bilato R., Hillairet J., Myra J., Kohno H., Poulos M. and Ono M. 2023 *Nucl. Fusion* **63** 026024

[25] Beers C.J., Green D.L., Lau C., Myra J.R., Rapp J., Younkin T.R. and Zinkle S.J. 2021 *Phys. Plasmas* **28** 093503

[26] Myra J.R. 2021 *J. Plasma Phys.* **87** 905870504

[27] Moritz J., Faudot E., Devaux S. and Heuraux S. 2016 *Phys. Plasmas* **23** 062509

[28] Jacquot J., Colas L., Clairet F., Goniche M., Heuraux S., Hillairet J., Lombard G. and Milanesio D. 2013 *Plasma Phys. Control. Fusion* **55** 115004

[29] Myra J.R. 2017 *Phys. Plasmas* **24** 072507
[30] COMSOL Multiphysics reference manual (available at: https://doc.comsol.com/6.1/docserver/)

[31] Lu L.F. et al 2023 *Nucl. Fusion* **63** 066023

[32] Lu L., Crombé K., Van Eester D., Colas L., Jacquot J. and Heuraux S. 2016 *Plasma Phys. Control. Fusion* **58** 055001

[33] Bobkov V. et al 2024 *PSI 26* (Marseille, France)

[34] D' Ippolito D.A. and Myra J.R. 2006 *Phys. Plasmas* **13** 102508

[35] Ochoukov R., Whyte D.G., Brunner D., D'Ippolito D.A., LaBombard B., Lipschultz B., Myra J.R., Terry J.L. and Wukitch S.J. 2014 *Plasma Phys. Control. Fusion* **58** 015004

[36] Hong R., Wukitch S.J., Lin Y., Terry J.L., Cziegler I., Reinke M.L. and Tynan G.R. 2017 *Plasma Phys. Control. Fusion* **59** 105008

[37] Colas L., Lu L.-F., Křivská A., Jacquot J., Hillairet J., Helou W., Goniche M., Heuraux S. and Faudot E. 2017 *Plasma Phys. Control. Fusion* **59** 025014

[38] Heuraux S., da Silva F., Ribeiro T., Despres B., Campos Pinto M., Jacquot J., Faudot E., Wengerowsky S., Colas L. and Lu L. 2015 *J. Plasma Phys.* **81** 435810503

[39] D' Ippolito D.A. et al 2008 *Phys. Plasmas* **15** 102501

[40] Kohno H., Myra J.R. and D'Ippolito D.A. 2015 *Phys. Plasmas* **22** 072504

[41] Urbanczyk G. et al 2024 *Nucl. Fusion* submitted

[42] Colas L. et al 2024 *Nucl. Mat. Energy* submitted