Validation of high-fidelity ion cyclotron range of frequencies antenna coupling simulations in full 3D geometry against experiments in the ASDEX Upgrade tokamak

We address the validation of finite elements ion cyclotron range of frequencies (ICRF) antenna coupling simulation against experiments performed in the ASDEX Upgrade tokamak. Measurements of the loading resistance in ICRF-heated, magnetically-perturbed 3D plasma discharges are compared against numerical predictions of the RAPLICASOL code. To this end, the 3D induction field and the 3D density profile are modeled by concatenating the PARVMEC, BMW and EMC3-EIRENE codes. The 3D density is input to RAPLICASOL, where full-wave simulations are performed on a finite element mesh retaining full 3D geometry in the ICRF antenna model and the plasma description. The results are further compared with RAPLICASOL simulations using a 1D density profile as measured at the outboard midplane in the same experiments. We find that simulations using a 1D density profile overestimate the change in loading resistance by a factor of ∼197 − 248%, while simulations using the full 3D density profile are in agreement with experiments within a factor ∼16%.


Introduction
Full-wave simulations of ion cyclotron range of frequencies (ICRF) antenna near fields in plasma are crucial to the understanding of power coupling to the core and mitigation of sheath further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. rectification effects. They are used to improve our knowledge about the underlying physical phenomena driving impurity sputtering and to design ICRF systems for future tokamaks, such as ITER and DEMO, that must outperform those in current machines. It is thus, essential, to validate these simulation tools against experiments in existing machines, which can point out flaws in terms of missing physics and made oversimplifications. One such simplification is that of assuming the plasma to be a slab, with density and induction field gradients only along the normal-to-the-antenna direction. While this approximation may work correctly under certain conditions, we will show in this paper that it does not correctly describe the antenna behavior in 3D plasmas. This affirmation will be seen true as long as the plasma perturbation scale-length is smaller than or of the order of the antenna poloidal and/or toroidal extent, and even when the deviation from axisymmetry is very small compared to the machine minor radius (~2 − 4% for ASDEX Upgrade).
In previous work [1,2], the loading resistance change of the four ICRF antennas installed in the ASDEX Upgrade tokamak was investigated in H-mode plasmas where an external, non-axisymmetric, magnetic perturbation (MP) field was applied and rigidly-rotated. In the explored lowcollisionality regime (ν * ei (ρ θ = 0.5) ∼ 0.07, ν * ei (ρ θ = 0.9) ∼ 0.54), the plasma response to the applied MP field consisted of increased radial transport across the separatrix, the so-called pump-out effect [3], and a 3D edge kink displacement. The latter, that we explore in this paper, is the result of amplification of the non-resonant components of the applied vacuum MP field by the plasma [4]. Such non-resonant equilibrium field leads to a pitch-aligned displacement of the flux surfaces from their original axisymmetric position [5]. For the sake of simplicity, we will assume the kink displacement part to be describable by ideal MHD theory. The effect of such flux-surface displacement on the coupling performance of ICRF waves can be better understood by introducing the concepts of R-cutoff and coupling region. The dispersion relation of plasma waves can be written in the cold plasma approximation as [6,7]: with n ∥ = ck ∥ /ω the parallel refractive index, where k ∥ is the component of the wave vector parallel to the local magnetic field, determined by the antenna current periodicity, R, L, S and P are elements of the cold plasma dielectric tensor: with ω the antenna angular frequency, Ω s = q s B/m s the cyclotron frequency and ω 2 p,s = n s q 2 s /ϵ 0 m s the plasma frequency of species s with density n s , charge q s and mass m s . The characteristic polynomial in equation (1a) yields as cutoff conditions for perpendicular propagation (n 2 ⊥ = 0): n 2 ∥ = R, the R-cutoff, n 2 ∥ = L, the L-cutoff, and P = 0, the cutoff for O-mode. In the ion cyclotron frequency range for a deuterium plasma, ω~Ω D , and for an on-axis magnetic field of B t ∼ 2.5 T, the R-cutoff occurs in the scrape-off layer (SOL) at n cutoff e ∼ 10 18 m −3 , the exact value depending on n 2 ∥ . The fast wave excited by the ICRF antennas in ASDEX Upgrade is thus evanescent in the volume bounded by n e < n cutoff e , where the antenna is embedded, and only propagates beyond it, where n 2 ⊥ > 0. The antenna near-fields thus depend on the optical conditions of this evanescent volume, also called coupling region, which likewise determine experimental quantities such as the reflection coefficient of transmission lines, the maximum voltage on them and the loading resistance [8][9][10][11][12][13].
The density perturbation produced by the kink mode excited when MPs are applied extends beyond the core region into the SOL. The R-cutoff layers of both the 2-strap antennas (k ∥ ∼ 7.5 m −1 ), and the 3-strap antennas [14,15] (k ∥ ∼ 11 m −1 in dipole phasing), as well as the evanescent region, acquire a homologous 3D deformation. We correlated the Rcutoff displacement locally registered by different low-field side (LFS) diagnostics (O-mode [16] and X-mode reflectometry [17,18], as well as the lithium beam [19]), to the change in loading resistance of the four antennas. An exponential fit of the type R L ∝ e −αd cutoff revealed for the MP discharges α values much smaller, α MP ∼ 7 m −1 , than those found when an axisymmetric plasma is radially moved (i.e. when the aspect ratio is scanned with changing major radius), α ∼ 18 m −1 . An example of this behavior is demonstrated in figure 1, where averaged absolute loading resistance measurements from the 2-strap antenna feeders are plotted against the absolute position of the R-cutoff, as determined by O-mode reflectometry. Suitable fits are included for the axisymmetric data (ASDEX Upgrade discharges #29 823, #30 629, #30 630 and #30 631 at t = 2 − 3.05 s) and the MP discharge (ASDEX Upgrade discharge #34 622 at t = 2.15 − 4.45 s). The fit uncertainty, denoted with a shadow, comprises a 99% confidence interval, assuming the data points normally distributed. The axisymmetric data was recorded with an upper plasma triangularity of δ u = 0.1, which best matches the poloidal antenna curvature. It is apparent that the loading resistance change is quite milder in the MP case, with very similar α values to the ones previously reported. The bifurcation in the loading resistance, when plotted against the limiter-R-cutoff distance, is expected since this is not the variable that dictates the optical properties of the coupling region. The optical thickness η =´k ⊥ dl is [9,10], but is of harder experimental access. Full-wave ICRF simulations with the RAPLICASOL code [20][21][22] using 1D density profiles, as measured by the lithium beam diagnostic, were seen to align with high-accuracy with the expected axisymmetric scaling, but failed to reproduce the change in loading resistance in the MP discharges.
In this paper, we perform RAPLICASOL simulations using instead a 3D density profile. To this end, we have first computed ASDEX Upgrade 3D magnetohydrodynamic (MHD) equilibria with MP coils using the PARVMEC code [24][25][26][27] for the MP discharges studied in [2] (#34 622, #34 634, #34 672, #34 673). This is done to validate against experiments the computed PARVMEC flux-surface kink displacements for a scan of poloidal field spectrum (∆φ UL ) with n = 2 toroidal periodicity in the MP coils. The computed plasma currents and flux-surface geometry for the ∆φ UL ∼ 0 • case are then used to extend the induction field domain to the SOL using the BMW code. The full 3D PARVMEC&BMW induction field is used in the EMC3-EIRENE code [28], where the 3D electron density can be reconstructed via comparison against lithium beam measurements. It is this 3D density that we use in RAPLICASOL, where a comparison of the obtained loading resistance against the experimental measurements is made. A workflow of the used codes in this paper and their relevant output quantities can be seen in figure 2. We find the new RAPLICASOL simulations using this 3D density profile to be in good agreement with the loading resistance measurements in MP discharges. This result highlights the need to treat the plasma density in 3D geometry even when small deviations from axisymmetry are present.
The paper is organized as follows: in section 2, we introduce the modeled discharge and the results from the PAR-VMEC code. A resolution scan over the PARVMEC input parameters is presented in appendix A. In section 3, we describe the extension of the magnetic induction field domain to the SOL using the BMW code. In section 4, the EMC3-EIRENE simulations and the comparison between modeled and measured 3D electron density are described. In section 5, RAPLICASOL simulations using the calculated 3D electron density are presented and compared to simulations using a 1D density profile and the experimental measurements of the loading resistance in MP discharges. In section 6, we draw the conclusions of the study.

3D ideal MHD equilibrium calculation
In this section, we describe the ideal MHD equilibrium reconstruction performed with the PARVMEC code. We compare the obtained equilibria against measurements of flux-surface displacement in the aforementioned set of MP ASDEX Upgrade discharges.

The PARVMEC code
The PARVMEC code is an MPI-capable version of the well known NEMEC code: VMEC (Variational Moments Equilibrium Code) plus NESTOR (NEumann Solver for TOroidal Regions). It seeks the minimum ideal MHD energy state of a given configuration, where the energy functional can be written as: Neumann boundary conditions for the normal component of the induction field and pressure continuity are included in order to compute free-boundary equilibria. The quality of convergence is often described through the ideal MHD force residual, which in straight field-line angle coordinates (s, θ * , ϕ) can be written as: With ψ ′ , χ ′ and p ′ the radial derivatives of the toroidal flux, poloidal flux and scalar kinetic pressure, √ g the Jacobian, H θ , H ϕ the covariant poloidal and toroidal induction field components and s the normalized toroidal flux, which serves as a flux-surface label. Here, ⟨·⟩ is the flux-surface average. The PARVMEC code was utilized with two objectives: first, as a way of numerically studying the plasma response to the applied MP field, and second, to provide a calculation of the 3D magnetic induction field for further simulations with the EMC3-EIRENE code. Since PARVMEC computes the resulting plasma equilibrium in the ideal MHD limit and Fourier representation, the code input is composed of the following elements: These profiles can be nowadays self-consistently integrated into PARVMEC directly from diagnostic data via the V3FIT code, which has also been recently coupled to SIESTA for the reconstruction of 3D plasma equilibria with islands [29]. Nevertheless, it is also possible to extract these quantities from a converged 2D equilibrium by another code. We have chosen the latter, as the former would require major effort integrating the ASDEX Upgrade diagnostics into V3FIT.

Selection of modeled discharge
First, a representative discharge and time point are chosen, from which the experimental data can be used as input for the axisymmetric reconstruction of the Grad-Shafranov CLISTE code [30]. Then, CLISTE's output can be used as the input for PARVMEC, and 3D effects can be simulated by further imposing the field from the MP coils. A discharge without MPs (the reference discharge in [2]), #34 632 at t = 5 s, is chosen as the basis for equilibrium reconstruction. This discharge acted as a 'reference' against which to compare the scan of MP field poloidal spectrum. The modeling results can be compared to measurements of the separatrix displacement since lithium beam, O-mode reflectometry and charge exchange recombination spectroscopy (CXRS) [31] measurements were available. The motivation to use a non-magnetically perturbed discharge instead of one with MPs relies on the simplification of diagnostic handling to obtain a more accurate CLISTE result. On the other hand, we are explicitly neglecting the effects of MPs on the plasma pedestal, such as the increased transport during pump-out. These effects will not arise in PARVMEC, since no transport equations are solved, and would have to be included indirectly by using MP-affected profiles as input. As it will be discussed at the end of this section, neglecting these effects might have an impact on the resulting simulations and how they compare to the experimental data. The input kinetic pressure and toroidal plasma current profiles for PARVMEC have been plotted in figure 3. Notice, that the pressure profile contains no flattening close to rational flux surfaces. This choice of p ′ ̸ = 0 at rationals will result in divergent Pfirsch-Schlüter currents, as expected from ideal MHD theory [32], increasing the solution sensitivity on the used radial grid. A scan over the number of flux surfaces, used Fourier harmonics and tolerance of the force residual in PARVMEC is presented in appendix A. The alternative choice of p ′ = 0 at rational flux surfaces was not explored in this work.
We now sketch the general steps that have been followed in performing the PARVMEC simulations. First, the CLISTE plasma boundary is fed into the DESCUR code [25], which computes its Fourier representation. Since the separatrix contains a singular point (X-Point) that cannot be represented by a finite number of Fourier harmonics, the toroidal flux is truncated just inside it, close enough that the loss of plasma volume is as small as possible. We have utilized ψ trunc ∼ 0.999 ψ sep CLISTE ∼ −3.15 Wb for the truncation. Next, the vacuum induction field is computed with the MAKEGRID code, which uses a Biot-Savart representation of the external conductors. Initially, we only compute the 2D induction field (no MPs), which will serve for the comparison between the CLISTE and PARVMEC free-boundary solutions in axisymmetry. Since CLISTE includes SOL currents, which are outside PARVMEC's computational domain, the first found solutions from the two codes will generally not coincide. In order to artificially include the effect of SOL currents in PARVMEC, and thus reduce the differences with respect to CLISTE, the following procedure is utilized: (i) The SOL currents are progressively decreased in the CLISTE solution, and the plasma boundary is newly evaluated each time.
(ii) The small differences in plasma boundary arising between the full CLISTE solution and the reduced SOL current one are compensated by adjusting (~1.5%) the currents in the external poloidal field coils. The procedure is repeated until the currents in the SOL have been reduced to a value which generates very small differences between the axisymmetric CLISTE and PARVMEC solutions. In praxis, once the CLISTE solution has been fixed, further small adjustments are made to the poloidal field coils taking the PARVMEC solution as a reference, as to improve the agreement even more. The final comparison of plasma boundaries, inner flux surfaces and qprofiles can be seen in figure 4. Here, we have used m θ = 28 poloidal harmonics, n s = 3001 flux surfaces and a tolerance of the force residual of f tol = 1 × 10 −15 for the PARVMEC simulation. The agreement between the two codes becomes excellent.
From this point on, an external 3D field from the MP coils is applied in the PARVMEC simulations, in order to evaluate the final 3D equilibrium. Since the MP coils are mounted on the passive stabilization loop (PSL), a bulk conductor meant to reduce the growth rate of the vertical instability, image currents develop which need to be taken into account. The PSL impact on the MP field poloidal spectrum and amplitude is evaluated by using the response function from finite elements calculations [33,34]. This effectively reduces the MP coil currents from the actuator value of I MP ∼ 5 kAt to I PSL MP ∼ 3.2 kAt. The plasma displacement, ⃗ ξ n , is computed from the final 3D equilibrium as the oscillation of a flux-surface position when intersected by a given plane. This is analogous to how diagnostics register oscillations of a given constant quantity pertaining to a given flux surface in the ideal MHD picture. Allegedly, this assumption discards the presence of magnetic islands, an approximation which so far has worked correctly when comparing VMEC predictions to experimental measurements of MP discharges [34].

Comparison to measured experimental displacements
A PARVMEC parameter scan over the number of Fourier harmonics, radial grid points and tolerance of the force residual is included in the appendix A. It was used to decide what were the best numerical input parameters, given the limitations imposed both by physics and numerics. In general, good convergence was found when the number of flux surfaces was greater than n s ≳ 501, the number of poloidal mode numbers was greater than m θ ≳ 30, the tolerance of the force residual was smaller than f tol ≲ 5 × 10 −14 , and the number of toroidal mode numbers was kept not larger than n ϕ = 4. In fact, further increasing n s and m θ does make the solution better converged (displacement and q-profile wise), but how much these parameters can be increased is limited by the available memory and wall-clock time in a given cluster.
Given these circumstances, a solution which involves a compromise between physical accuracy and numerical resources was chosen, i.e. n s = 1251, m θ = 32, n ϕ = 4 and f tol = 5 × 10 −14 . The vacuum field cylindrical grid was kept at {n r = 256, n z = 512, n ϕ = 64} points with 2 toroidal field periods. A scan over the MP field poloidal spectrum, ∆φ UL , is made with this accuracy and the PSL-attenuated coil currents, I PSL MP ∼ 3.2 kAt. Since the MP field rigid rotation is performed with ν = 3 Hz and n = 2 toroidal symmetry, we compared the ν = 3 Hz fundamental sine amplitude of an ordinary least squares (OLS) fit to the separatrix position (d sep ) measured in the experiments reported in [2] of the type: and the n = 2 component of the PARVMEC simulations. Here, | ⃗ ξ n (ν j )| is the sine amplitude of the OLS fit for the harmonic with frequency ν j and phase ϕ j , α k are polynomial coefficients and t is the time variable. Three harmonics were used for the sine fitting with a cubic polynomial for baseline correction.
The result of such comparison can be seen in figure 5(a). It is apparent that the numerically computed displacements fall short from the experimental ones. Several factors may contribute to this disagreement, mainly: (a) Uncertainties in the pressure and current profiles of the underlying CLISTE equilibrium, for instance, related to core ICRF fast ions. (b) Lack of pump-out in the pedestal profiles. The CLISTE solution belongs to discharge #34 632, the one without MPs. Therefore, no pump-out effect is present in the kinetic profiles. Since PARVMEC preserves flux-surface quantities, this effect is not reproducible within the simulation when the MP field is imposed, and needs to be included at the input profile level. (c) Insufficient numerical accuracy in the PARVMEC solution. As stated previously, increasing further m θ and n s does increase the n = 2 component of the OMP displacement non-negligibly. (d) Uncertainties in the PSL model. The finite elements calculations in [33] are computed for a single coil-plasma distance, making them a first-order approximation. (e) Ideal MHD approximation. If the real magnetic equilibrium contains islands which width is determined by resistive MHD physics, this cannot be captured with the used model.
For these reasons, a scaling factor is used for the currents of the MP coils. The mean of the experimental displacements is computed for each ∆φ UL . The ∆φ UL = −2.2 • is used to derive the factor, such that: with ξ exp n (ν = 3 Hz) ∼ 8.98 mm and ξ PARVMEC n (n = 2) ∼ 4 mm. Each MP coil current is scaled by this factor, yielding good agreement between PARVMEC and the experimental measurements for all the discharges. The addition of said factor is further motivated in this case, because PARVMEC is used as a reconstruction tool, rather than as a predictive model. We compare in figure 5(b) the CXRS separatrix position for the ∆φ UL = −2.2 • case overlaid for several MP periods and the PARVMEC last closed flux surface (LCFS). Both are found to be in very close agreement, aside from a relative radial shift of ∆R ∼ 1 cm. An OLS fit reveals that the n = 2 component phase shift between the scaled PARVMEC equilibrium and the CXRS time traces is of the order of ∆ϕ~0.6 • (∼ 0.01 rad). In principle, this negligible phase shift would reinforce the validity of the used PSL model. We note here the previous statement that the deviation from axisymmetry is of the order of 1 cm in a ∼ 50 cm (ASDEX Upgrade minor radius), that is 1/50~2%, or 4% if taken peak to peak. Despite this low deviation, we will see in section 5 that an axisymmetric approximation to the density used for antenna near field simulations fails to reproduce the experimental loading resistance measurements.

Full-domain magnetic induction field computed with the BMW code
The BMW (Biot-Savart Magnetic VMEC Vector potential) code calculates the magnetic vector potential arising from the plasma current, ⃗ ȷ p , by a volume integral: and then obtains the 'plasma' magnetic induction field as The total induction field is thus defined as the sum of plasma and vacuum fields, the latter produced by the external conductors: BMW is directly coupled to PARVMEC by using the computed plasma currents and flux-surface geometry for the volume integral. This procedure is applied to the core and SOL regions thus allowing to extend the resulting equilibrium outside the LCFS, which is a fundamental step for SOL transport studies. Due to this curl approach, BMW produces inherently divergence-free ⃗ Bfields, since ⃗ ∇ · ⃗ ∇ × ⃗ A p = 0 and ⃗ B v is divergence-free from the Biot-Savart calculation. This feature allows its usage with transport codes such as EMC3-EIRENE, which depend on a pre-calculation of the background magnetic induction field that must be continuous across the separatrix, in order to avoid non-physical behavior. In this section, we detail the main topological characteristics of the magnetic induction field chosen for the forthcoming steps in our simulation workflow. The scaled-currents, ∆φ UL = −2.2 • PARVMEC simulation previously described was extended with the BMW method to obtain the full-domain induction field including the SOL. The cylindrical grid where the magnetic vector potential is computed was kept at the same resolution used for PARVMEC, i.e. {n r = 256, n z = 512, n ϕ = 64} with 2 field periods. The resulting connection length 7 at the OMP was compared against chosen PARVMEC flux surfaces as seen in figure 6. The agreement is impoverished at the very edge, where islands form densely due to the elevated shear. This is due to the expected divergent Pfirsch-Schlüter currents at rational flux surfaces with p ′ ̸ = 0 and the stringent numerical requirements that this situation poses to the PARVMEC solution (see appendix for a detailed discussion). Even if the resulting PAR-VMEC induction field is nested, this is only up to the chosen simulation accuracy. When the computed plasma currents are employed in the aforementioned method, these are not able to completely shield the external vacuum field B ψ component, yielding islands. However, if we compare a non-rational flux surface (here chosen at ρ θ = 0.87), both the surface and the high-connection length region overlap satisfactorily. This proves the reliability of the method to preserve the magnetic topology from PARVMEC where no rational flux surfaces exist, while in the domains where they do exist, it is generally not possible. As an example, we compare the contour line for L c = 100 m with the PARVMEC LCFS, which quickly exposes that these boundaries are topologically different. This effect can be further verified with a Poincaré plot of the BMW magnetic induction field mapped to the {ρ θ , θ} plane. Field lines were followed with the GOURDON code [35]. The poloidal flux is computed in axisymmetric approximation from the BMW toroidal component of the magnetic vector potential: where χ sep and χ axis are taken to be the interpolated poloidal flux values at the separatrix and magnetic axis of 7 The length along a magnetic field line between the point it is started from and its intersection with a wall structure. the underlying PARVMEC solution. The presence of islands becomes apparent for all rational flux surfaces with m, n = 2 topology as can be seen in figure 7. Further islands can also be discerned, such as the m = 7, n = 4, m = 9, n = 4, etc. The position of the rational flux surfaces from PARVMEC has been included as red vertical lines, denoted by their rational number. The shift between the actual island and the PARVMEC rational flux surface position is expected, due to the implicit ι-jump that is produced for a given resolution in order to preserve nestedness [36] and the approximation in equations (8) and (9). Despite the appearance of numerical islands in the equilibrium core, we will see in the next sections that this does not influence our intended study, which happens mostly in the SOL.

EMC3-EIRENE kinetic profile reconstruction
The EMC3-EIRENE transport calculations can be summarized in two steps. One of technical nature and the second of more physical relevance. The first one is the construction of a computational mesh given that the BMW solution contains closed flux surfaces, islands and stochastic regions. The second one is the appropriate choice of perpendicular transport coefficients for particle (D ⊥ ), electron heat (χ e ) and ion heat (χ i ) fluxes. We describe both steps in what follows.

Grid geometry
In this step, a 3D grid is constructed taking as a reference the PARVMEC&BMW magnetic induction field solution described in the previous section. First, a 2D grid is built in each toroidal plane. This grid is usually meant to be aligned with the separatrix as best as possible, in order to provide sufficient radial resolution in the pedestal region. In our field solution, islands overlap at the very edge, making this step not trivial. The inner boundary of the simulation domain in a given 2D grid is made to align with a non-rational flux surface, which provides a closed surface avoiding particle losses. Afterwards, magnetic field lines are followed from the nodes of each 2D grid until the next one, thus building the final 3D grid. Flux tubes from which field lines can be bi-linearly interpolated at runtime are constructed in this step. This procedure saves significant computational resources, avoiding a real-time field line tracing for each Monte Carlo particle. A simulation with no sources and no fluxes through the boundary surfaces is performed to check the quality of the grid, the so-called Monte Carlo test. If the grid quality is optimal, a uniform particle probability density is expected. Furthermore, flux conservation within each grid generated flux tube is checked to further diagnose the quality of the grid. The 2D grid geometry with an overlaid Poincaré plot, as well as the result from the Monte Carlo test and flux conservation diagnostic can be seen in figure 8. The probability distribution function displays random fluctuations with an average value of ⟨ f density ⟩ = 1.01 ± 0.07, where the error represents 3σ. The fact that no particular geometric feature can be associated with such fluctuations, and their small standard deviation, is proof of the good divergence-free nature of the magnetic field. Likewise, deviations to flux conservation amount only to max(δ flux ) ∼ 0.3% (except very close to the MP coils), in good agreement with previous studies of EMC3-EIRENE grid formation from a PARVMEC&BMWC equilibrium in the W7-X stellarator [37]. Given the good quality of the grid diagnostics, we proceed to use it in the calculation of plasma kinetic profiles.

Selection of transport coefficients
An optimum set of transport coefficients is that which best reproduces the experimental results. To verify the agreement, the lithium beam electron density data has been compared against the EMC3-EIRENE simulations. The computed n e is mapped to the lithium beam LOS and compared against the experimentally reconstructed n e profiles. Usually, it is desirable to compare against more than one diagnostic, for instance, thermionic divertor currents [38], electron cyclotron emission for T rad profile [39], etc can be included. In this case, only the lithium beam was used, and additional ones will be left as future work. In order to quantify the degree of agreement between the modeled and the experimental data, we define the following metric: which shall be minimized. Here δ min/max is the metric at ϕ = ϕ min/max , i.e. where the kink displacement in the lithium beam LOS is minimum/maximum and x i are the lithium beam channel positions in the LOS, such that 10 16 m −3 < n exp e < 2.5 × 10 19 m −3 . This density range is chosen because it covers completely the density region used in the simulation domain of RAPLICASOL, as it will be described in the next section. The reconstruction of kinetic profiles was performed as follows: A 1D cross-field particle diffusion coefficient profile, D ⊥ , was chosen as a function of ρ θ . This profile includes a reduced transport region inside the LCFS, which is representative of the core high-confinement properties, a linear increase in the near SOL, where particles are no longer well confined, and a constant layer for the far SOL. This 1D D ⊥ profile is then mapped to the EMC3-EIRENE 3D grid by using the ρ θ coordinate from BMW (equation (9)). After testing different D ⊥ profiles, it turned out that a step in the far SOL needed to be introduced. The lithium beam is toroidally close to the limiter of a 2-strap antenna, which likely affects the density profile in the far SOL. However, in EMC3-EIRENE, even though realistic limiters are included in the simulation, these are fixed with respect to the MP field rotation phase, i.e. we are simulating a single time point (and then we will just rotate the electron density) instead of performing a new EMC3-EIRENE simulation for each MP rotation phase. Therefore, the modeled electron density is not affected at the maximum and minimum LOS kink displacements by this limiter, unlike the lithium beam diagnostic is in the experiment. This is, however, not of great importance since we set a vacuum layer in RAPLICASOL for densities lower than 2 × 10 17 m −3 < n e anyway. As a boundary condition, the electron density was set to n e = 2 × 10 19 m −3 just inside the separatrix. The resulting n e , T e , D ⊥ and neutral particle deuterium density can be seen in figure 9 for a poloidal cross-section at ϕ~0. It is readily seen that the density and temperature profiles are 3D perturbed. In particular, the characteristic lobes from the application of MPs can be more clearly seen in the X-Point region of the temperature profile. The comparison between the modeled and experimental electron density can be seen in figure 10. The resulting metric for this case is δ tot ∼ 2.1 × 10 18 m −3 . The EMC3-EIRENE n e profile is taken at the lithium beam LOS for ϕ min and ϕ max and compared against the experimental profiles, which have been obtained by the phase-lock average of n e over all MP periods. Very good agreement is achieved for 10 18 m −3 < n e < 1.5 × 10 19 m −3 . The profile in the far SOL (n e < 10 18 m −3 ) decays slower than the measurements, due to the lack of the limiter effect. The profile inside the LCFS n e ⪆ 1.5 × 10 19 m −3 flattens out rapidly, which points to the influence of the edge stochastic region in the PARVMEC&BMW solution. Likewise, the same effect can be observed when the 2D n e in the lithium beam LOS is compared, as seen in figure 11. The agreement in the SOL is very good, but the simulated core density is lower due to the increased transport in the stochastic region. We also observe the ∼ 1 cm radial shift of the separatrix density reported in figure 5(b). Overall, the result is very successful, since: (a) An EMC3-EIRENE grid has been built from the PAR-VMEC&BMW induction field solution for a tokamak for the first time, that does not suffer from the violation of ⃗ ∇ · ⃗ B ̸ = 0 across the LCFS. (b) The reconstructed kinetic profiles in the SOL display the characteristic topology expected from the application of a 3D MP field. Furthermore, good agreement is found when compared to the experimental lithium beam profiles for the discharge we intended to model. (c) The R-cutoff displacement of the fast wave matches well that of the experiment, which validates the usage of this 3D profile for the forthcoming RAPLICASOL simulations.

3D RAPLICASOL ICRF coupling simulations
In this section, we perform RAPLICASOL ICRF coupling simulations for the 2-strap antenna with the computed EMC3-EIRENE 3D density profile. The magnetic induction field is treated in 1D approximation, due to performance limitations in the absorbing boundaries (perfectly matched layers, PMLs) with 3D magnetic fields 8 . We perform a set of 7 simulations by toroidally rigidly rotating the density profile. The toroidal center of the ICRF antenna was thus placed at The intention is to reproduce the same coupling variations found in experiments when the MP field is rotated. In figure 12, we show the antenna center positions along the 3D density R-cutoff at the OMP, and the 1D radial n e profiles taken at those same locations. A least squares fit to the R-cutoff position is performed, like previously done for the separatrix (see equation (3)). The plasma domain used in RAPLICASOL is comprised between n min e = 2 × 10 17 m −3 ≤ n e ≤ n max e = 1.4 × 10 19 m −3 . Here, n min e is the minimum density value allowed, which prevents the appearance of the lower hybrid resonance. Any density below this value is assumed n e = n min e up to the vacuum layer where the antenna is embedded. n max e is the maximum density allowed, which prevents gradients in the plasma-PML interface. Any density above this value is assumed n e = n max e . This value is somewhat smaller than that used in the simulations using a 1D density profile reported in [2], in order to avoid the artificial density-well inside the LCFS, which arises at the edge stochastic region. For the new simulations, we use uniquely a curved antenna model of the 2-strap antennas based on experimental measurements of their limiter positions and a CAD design of the antenna. The antenna model, as well as the PMLs, can be seen in figure 13. The RAPLICASOL antenna loading resistance change is evaluated by applying the same least squares fit to the obtained per-port loading resistances, and then taking the average change from both ports: Where ⟨R L ⟩ is the average loading resistance of the 7 simulations, and ∆R L is the n = 2 amplitude of the least squares fit. The resulting per-port loading resistances can be seen in figure 14(a),(b). Their average is compared to RAPLICASOL simulations using 1D density profiles, as measured by the lithium beam diagnostic, and the experimental loading resistance time traces for the modeled discharge, #34 622 at t = 2.15 − 4.45 s in figure 14(c). Since the EMC3-EIRENE 3D density reconstruction is validated against the lithium beam diagnostic, only the 2-strap antenna toroidally closest to it, the so-called ICRH3~16 • apart, is displayed. The second 2-strap antenna is~164 • apart from the lithium beam, making the reconstructed density less accurate in that region, because of toroidally localized plasma asymmetries. The time coordinate in experiments is mapped to the toroidal angle in the simulations via the formula: where ν MP = 3 Hz, n MP = 2 and t 0 = 3.5 s is the time point used to read in the MP coil currents from the experiment, which defines the toroidal phase of the MP field and that of the density locked to it. Due to uncertainties in the absolute    radial positions of the plasma and ICRF antenna, the loading resistance from RAPLICASOL has been multiplied by a factor to equalize the means. The experimental average is taken over the two feeders of ICRH3 and all the MP cycles. We observe that the average loading resistance is very well reproduced by the RAPLICASOL simulations employing the EMC3-EIRENE 3D density profile, which was the objective of the study. RAPLICASOL simulations employing a 1D density profile, on the other hand, greatly overestimate this average. Note, however, that the α exponent from simulations using 1D density profiles, α ∼ 18.5 m −1 , is in excellent agreement with the axisymmetric radial scan displayed in figure 1, as would be expected (the poloidal plasma shape was matched to that of the antenna limiter, thus minimizing poloidal and toroidal inhomogeneities).
Having obtained the n = 2 change of both the cutoff displacement and the loading resistance, we can further compare their relation to an experimental scaling elaborated with all the MP phasings displayed in figure 5(a), and a RAP-LICASOL scaling using the aforementioned lithium beam 1D density profiles, both reported in [2]. A diagram is presented in figure 15. The experimental scaling is derived by using the averaged α values obtained from the combined 2-strap antennas, ICRH1&ICRH3 data sets, resulting in α 2−strap = 7.22 m −1 (quite similar to the value displayed in figure 1). We also include the resulting scaling when using the data from ICRH1 alone, α ICRH1 = 8.64 m −1 , and that of ICRH3 alone, α ICRH3 = 5.73 m −1 . We have also included an uncertainty region comprised between the smallest measured scaling α(ICRH3) = 5.30 m −1 using the R-cutoff obtained from O-mode reflectometry, and the largest α(ICRH1) = 9.04 m −1 , obtained from the lithium beam, which covers both individual 2-strap antenna data sets together. The difference, albeit small, between the ICRH1 and ICRH3 scalings highlights the importance of locally resolving the electron density in the antenna vicinity. The experimental loading resistance median over all ICRH1 and ICRH3 feeders is further included as a scatter, as a function of the R-cutoff displacement measured by lithium beam, O-mode and X-mode (X4) reflectometry. It is readily verified that the RAPLICASOL prediction with a 3D density profile is more accurate than that previously obtained with 1D density profiles. This corroborates the nonaxisymmetric impact on the loading resistance, which cannot be well predicted in 1D approximation. The upper feeder of the computed RAPLICASOL loading resistance change with the 3D density profile falls within the uncertainty region, but the average value is just short. As a metric for comparison, we list the OLS α regression coefficients from experiments, RAPLICASOL simulations with a 1D density profile and those with the EMC3-EIRENE 3D density profile in table 1. These also include RAPLICASOL simulations using a flat antenna model (not discussed in this paper). In the 1D density profile simulations, ∆d cutoff corresponds to the R-cutoff displacements of the used 1D lithium beam density profiles, either directly measured for different MP phasings, or radially shifting one of the profiles, denoted as 'radial scan'. We also list the α exponent relative error between simulations and experiments, where the scaling of ICRH3 has been taken as the comparison point since the density reconstruction is most relevant for this antenna. When comparing the RAPLICASOL data sets with that of ICRH3, the relative error is reduced to a~16% for the 3D density profile, while 1D density profiles result in~197%-248%.

Conclusions
We have performed 3D ICRF antenna coupling simulations based on the computational workflow presented in figure 2. This addresses (i) the calculation of an ideal, non-linear, freeboundary MHD equilibrium for the core region of ASDEX Upgrade under applied MPs with the PARVMEC code; (ii) the extension of the magnetic induction field domain to the SOL via a volume integral of the core plasma currents with the BMW code; (iii) the preparation of a computational grid for the EMC3-EIRENE code, which served to verify the optimum quality of the magnetic field calculation; (iv) the reconstruction of 3D plasma kinetic profiles via comparison to experimental data and (v) the calculation of the loading resistance change, as the EMC3-EIRENE 3D density profile is rigidly rotated, with full-wave simulations using the RAPLICASOL code.
A converged 2D equilibrium from the CLISTE code was used as a starting point for the 3D PARVMEC simulations. A sensitivity scan was performed, which allowed identifying the main numerical parameters to optimize. We obtained that the used poloidal mode numbers, m θ , and radial grid points, n s , should be increased as much as possible to properly resolve the edge q-profile and the Pfirsch-Schlüter currents at rational flux surfaces. Meanwhile, the toroidal mode numbers, n ϕ , should not be abused, especially when they cannot be well represented at the plasma edge in simulations with elevated qprofile, like tokamak ones. The tolerance of the force residual, f tol , was found to play a role only above a certain threshold, i.e. f tol ∼ 5 × 10 −14 for this case. With a solution based on a compromise between numerical resources and physics fidelity, the computed plasma displacements fell a factor~2 short with respect to those measured in experiments. The main reasons are: (i) uncertainties in the input profiles, (ii) lack of a pumpout effect in the 3D simulation, (iii) limited numerical resolution and (iv) limitations of the ideal MHD approximation. In order to overcome this, the currents in the MP coils were multiplied by a factor f MP = 2.25 to match numerical and experimental displacements.
Through the extension of the magnetic field domain to the SOL with the BMW code, we verify that the flux-surface displacement from the PARVMEC solution is preserved for domains with no rational flux surfaces. This allows reproducing experimentally-relevant particle transport and the reconstruction of kinetic profiles with the EMC3-EIRENE code. Plasma domains with rational flux surfaces, however, lead to magnetic islands due to the imperfect shielding arising from the limited numerical resolution in PARVMEC. This is especially critical at the plasma edge, where a high magnetic shear exists.
An EMC3-EIRENE grid could be constructed, which verified the good quality of the induction field solution in terms of ⃗ ∇ · ⃗ B = 0. This grid is used to compute plasma kinetic profiles by adjusting the cross-field transport coefficients, D ⊥ in particular. The EMC3-EIRENE modeled electron density is compared against measurements of the lithium beam diagnostic. Their agreement is evaluated through a metric, δ tot , used to determine the optimal choice of transport coefficients. A satisfactory 3D n e profile could be reconstructed within the most relevant portion of the ICRF coupling region. The agreement impoverished close to the limiter region, where the lithium beam is affected by toroidally localized vessel structures that are not taken into account when one toroidally rotates a single EMC3-EIRENE simulation. Beyond n e ⪆ 1.5 × 10 19 m −3 , the results also depart from experiments, due to the stochastic edge in the magnetic field solution.
Full-wave simulations with the RAPLICASOL code have been performed using as input the EMC3-EIRENE 3D n e profile. The curved 2-strap ICRF antenna model was exclusively used. The density profile was toroidally rotated in several increments, and the loading resistance change was related to the R-cutoff movement. The resulting scaling is in better agreement with experiments within the considered experimental uncertainty. Unlike simulations employing a 1D n e profile, where an overestimation by a factor 197 − 248% is present, simulations employing the 3D n e reconstruction approximate the experimental results within a factor~16%.
Overall, despite the challenging experimental conditions, (i.e. characterization of plasma displacements and ICRF coupling changes in 3D geometry, which prevents the usage of common tools that assume axisymmetry), and large numerical challenges, i.e. development of a state-of-the-art 3D workflow, the found agreement between experiments and simulations is remarkable. These results settle the necessity to treat ICRF 3D coupling studies in full 3D geometry, as reduction to 1D results too crude to give predictive capabilities. The developed tools are validated, modular and flexible. These can be adapted to different tokamak and stellarator ICRF systems for their better future optimization. achieve. Furthermore, the vacuum field resolution (in a cylindrical grid, {n r , n z , n ϕ }, i.e. number of radial, vertical and toroidal points) can be investigated. It is also possible to optimize other computational parameters, such as the number of calls to the NESTOR solver, but these mostly influence the numerical convergence and are of less physical relevance to the result. A 3D run with ∆φ UL ∼ 0 • and PSL attenuated currents in the MP coils was chosen for the forthcoming scans, as it presents large enough LCFS displacements to more easily appreciate the effect of scanning a given parameter.

Spectral resolution:
In order to resolve the basic helical geometry of the plasma boundary with q LCFS ∼ 7.5 and an n = 2 MP toroidal field symmetry, we need at least m θ = 15 poloidal mode numbers in a straight field-line angle representation. Higher-order toroidal harmonics would consequentially require a larger number of poloidal harmonics in order to be represented at q LCFS . Since the number of MP coils is finite, the perturbation field does indeed contain such harmonics. In first approximation, the coils can be modeled as a square function in toroidal direction, which discrete Fourier series contains odd harmonics, n = 1, 3, 5.. (the equivalent continuous transform is sinc(πϕ)). In praxis, however, the number of m θ that can be used is limited due to the memory required to compute their Fourier coefficients. Two resolutions of n ϕ = 4 (n = −8, −6, ..., 6, 8 due to the used 2-fold toroidal symmetry) and n ϕ = 8 (n = −16, −14, ..., 14,16) were utilized in PAR-VMEC to study the influence of these high-order toroidal harmonics on the resulting equilibrium. The number of m θ was scanned, while n s , f tol , and the vacuum field resolution were kept constant. We used two numerical diagnostics to evaluate the quality of the equilibrium convergence: the full plasma displacements and n = 2 component at the outboard midplane, the plasma top, the inner midplane and the 'X-Point', as well as the q-profile at the LCFS. In ideal MHD, the q-profile is a topological invariant [40], and therefore, equilibria with an initial q-profile can only be converted through energy minimization into equilibria with the same q-profile. Would the q-profile not coincide, it would indicate that the used resolution is insufficient. The result of the comparison can be seen in figure A1.
Both cases were converged with n s = 1001 flux surfaces; the equilibria with n ϕ = 4 was converged down to f tol = 5 × 10 −14 , while n ϕ = 8 down to f tol = 5 × 10 −11 (smaller tolerance was often not reached for this case). Both cases display increasing (and convergent) displacements at all locations as the number of used m θ increases. However, two different behaviors are observed for q LCFS . In the n ϕ = 8 case, q LCFS ∼ 6.85 < q CLISTE LCFS ∼ 7.5 ∀m θ , while for n ϕ = 4, q LCFS → q CLISTE LCFS as m θ becomes larger. This indicates that the n ϕ = 8 resolution is not numerically stable, as the q-profile is not preserved in the energy minimization, and the increasing poloidal resolution does not solve the problem. We should keep in mind, that at this resolution the pitch-aligned higher toroidal mode numbers are only well represented at the most inner flux surfaces, i.e. for q = 6, n = 6 already precises m = 36, which is not reached in the scan. More information can be obtained after inspection of the edge flux surfaces, which display a clear magnitude increase in the high-order toroidal harmonic amplitudes for this case, displayed in figure A2. In the n ϕ = 4 case, increasing m θ leads to a better-resolved equilibrium, while the contribution of higher-order toroidal mode numbers remains small. In the n ϕ = 8 case, however, an oscillation of the plasma boundary becomes apparent at the LCFS. Increasing further m θ > 32 only increases the deviation from the expected n = 2 behavior. A similar observation has already been made in the literature. In [41], VMEC 9 was used to study n = 3 perturbed DIII-D equilibria. The EFIT axisymmetric code provided the necessary input in the same fashion as CLISTE is used in our study. It was found that when the EFIT equilibrium was truncated close to a rational flux surface, the 3D equilibria constructed from it would display poor convergence with respect to f tol and similar boundary oscillations of the flux-surface averaged toroidal current density at the LCFS. The suggested explanation is the ideal MHD instability of equilibria with near-separatrix rational flux surfaces. The proposed solution was to truncate the equilibria far enough inside the separatrix of EFIT, such that a large enough gap between main n-number rational flux surfaces and the boundary exists, but not deep enough that a significant fraction of the bootstrap current would be neglected. This method, however, presents several problems: (i) Rational numbers are dense in real space, thus deeming the term 'gap between rationals' subjective. (ii) This truncation often requires removing a significant fraction of the q-profile, since close to the separatrix rational flux surfaces become densely packed. For instance, in the aforementioned publication q trunc ∼ 5 as opposed to the EFIT value close to the LCFS of q EFIT LCFS ∼ 7.5. The resulting equilibria can therefore not represent correctly the poloidal structure of the edge displacements (i.e. with q trunc ∼ 5 and n = 3 implies m trunc ∼ 15), while the real poloidal mode number would be closer to the original EFIT edge value of q EFIT LCFS ∼ 7.5 → m EFIT LCFS ∼ 22, if we assume the kink displacement to be pitch-aligned at the LCFS. Since our aim is to correctly reconstruct the poloidal and toroidal mode structure in order to accurately describe the density perturbation in the LFS, this method could make our calculations less experimentally relevant.
The chosen approach in our case was to keep as much as possible of the toroidal flux by avoiding truncation deep inside the separatrix of CLISTE. This means limiting the number of used toroidal harmonics to n ϕ = 4, which results in a numerically stable and convergent result. Furthermore, this resolution is enough to resolve the two main discrete harmonics of the MP coils in 'square approximation', i.e. n = 1, 3 (n = 2, 6 for two field periods).
The first term on the rhs. represents the homogeneous solution, with f (ψ) an arbitrary function, and the second the inhomogeneous one. By plugging equation (12) into (), we obtain an expression for the parallel current density: The first term on the rhs. corresponds to the force-free δ current aforementioned. The second one, is the Pfirsch-Schlüter current. This last current is proportional to the pressure gradient, and can be seen to become singular at ι = n/m rational surfaces. In order to provide appropriate ideal MHD force balance with p ′ ̸ = 0, these currents need to be numerically resolved, otherwise the final equilibrium will fail to shield resonant external perturbations. The development of these currents in PARVMEC has been studied in detail in a simplified screw pinch geometry [36], and recently, in general toroidal geometry [46], leading to the conclusion that these are asymptotically resolved as resolution increases (in this case, the number of radial grid points).
A scan over the number of radial grid points was performed with the case n ϕ = 4, m θ = 34 for n s = 501, 1001, 1501, 2001. All simulations were converged down to f tol = 5 × 10 −14 . The displacements at the OMP, plasma top, HFS midplane and near the X-Point are shown in figure A3. We see that the change in plasma displacements is nearly independent of the number of flux surfaces, provided enough are used, i.e. n s ⪆ 501. It was shown in section 3, however, that even high radial grid densities of the order of n s = 1251 cannot provide enough resolution as to completely shield experimentally-relevant resonant external perturbations.
Tolerance of the force residual: This scan intends to estimate a reasonable convergence criterion in terms of the force residual described in equation (2). The main limiting factor to minimize this residual is the wall-clock time of the cluster (i.e. the maximum time a simulation is allowed to run), and the gradient of the plasma thermal energy minimization, which generally, decreases considerably the smaller the force residual becomes. The case with n ϕ = 4, m θ = 34, n s = 1001 was converged down to f tol = 1 × 10 −10 , 1 × 10 −12 , 5 × 10 −14 , 1 × 10 −16 . The full and n = 2 component of the plasma displacement at the OMP, plasma top, inner midplane and X-Point was used as a diagnostic tool for convergence, and the result is displayed in figure A4.
The case with f tol = 5 × 10 −14 is practically indistinguishable from the one with f tol = 1 × 10 −16 . Therefore, the former limit has been used in the forthcoming simulations. Another interesting result is the spatial distribution of the force residual, which can be locally computed following equation (2). This is displayed in figure A5, where the toroidally-averaged local force residual has been normalized to the local pressure gradient. It is seen that, even though the innermost flux surfaces tend to relax asymptotically the smaller f tol becomes (f local tol ∼ 10 −8 | ⃗ ∇p|), the outermost ones do not (f local tol ∼ 10 −3 | ⃗ ∇p| − 10 0 | ⃗ ∇p| ∀f tol ). This is clear from a numerical and physical perspective. It is in this region where most rational flux surfaces exist in the smallest volume, and therefore, where most numerical accuracy is required in the radial grid. In turn, this also means that for a given accuracy in the radial grid, number of Fourier harmonics and poloidal/toroidal points, further minimization of f tol will not significantly improve the solution at the edge, where more numerical accuracy may be required instead.