Iterative addition of parallel non-local effects to full wave ICRF finite element models in axisymmetric tokamak plasmas

The current response of a hot magnetized plasma to a radio-frequency wave is non-local, turning the electromagnetic wave equation into an integro-differential equation. Non-local physics gives rise to wave physics and absorption processes not observed in local media. Furthermore, non-local physics alters wave propagation and absorption properties of the plasma. In this work, an iterative method that accounts for parallel non-local effects in 2D axisymmetric tokamak plasmas is developed, implemented, and verified. The iterative method is based on the finite element method and Fourier decomposition, with the advantage that this numerical scheme can describe non-local effects while using a high-fidelity antenna and wall representation, as well as limiting memory usage. The proposed method is implemented in the existing full wave solver FEMIC and applied to a minority heating scenario in ITER to quantify how parallel non-local physics affect wave propagation and dissipation in the ion cyclotron range of frequencies (ICRF). The effects are then compared to a reduced local plane wave model, both verifying the physics implemented in the model, as well as estimating how well a local plane wave approximation performs in scenarios with high single pass damping. Finally, the new version of FEMIC is benchmarked against the ICRF code TORIC.


Introduction
Ion cyclotron resonance heating (ICRH) is a robust and versatile method for heating fusion plasmas that is commonly used in present fusion experiment and planned to be installed in both ITER and SPARC [1,2].In large scale devices it is the only heating scheme available that can provide dominant ion heating, making it an attractive system for heating the plasma to fusion relevant conditions in a reactor.In addition, ICRH can be used to reduce turbulent transport, control impurities, and to provide non-inductive current drive [3][4][5][6][7].
The most common use of ICRH is to launch a fast magnetosonic wave that is absorbed by ions at fundamental or harmonic cyclotron resonances, or by the electrons via a combination of transit time magnetic pumping (TTMP) and electron Landau damping (ELD) [8].In addition, the wave may be mode converted into ion Bernstein waves (IBW) and ion cyclotron waves (ICW) [9][10][11].Although, if the resonance and ion-ion hybrid layers overlap, which they do for standard minority heating scenarios in large and most medium-sized tokamaks, mode conversion is negligible [12].In this work we will therefore focus on numerical modelling of the fast magnetosonic wave and not consider mode conversion.
At fusion relevant temperatures and densities, the mean free path is longer than the gradient scale lengths.The dielectrically induced current is therefore non-local, meaning that the velocity of a particle in a given point is dependent on the cumulative acceleration in all previous positions.This turns the electromagnetic wave equation into an integrodifferential equation, which is not natural to solve using standard numerical methods, such as the finite element method (FEM).Furthermore, the wavelengths of the fast magnetosonic wave are non-negligible compared to the size of the machine, indicating that the wave equation is most appropriately solved using a full wave solver.
The non-local nature of the plasma exhibits different properties perpendicular and parallel to the magnetic field lines.In the perpendicular direction, the non-locality arises due to finite Larmor radius (FLR) effects.The particle makes excursions from the magnetic field line and will be subject to different electric field strengths along its gyro-orbit.In fast wave heating scenarios, the ion Larmor radius is usually small compared to the wavelength and is therefore often used as a smallness parameter for expansion of the perpendicular plasma response tensor.Using such an expansion, the perpendicular non-locality can be modelled using for example a FEM description [13][14][15][16].Perpendicular non-locality is responsible for effects such as higher harmonic cyclotron damping, TTMP, and mode conversion to the IBW [17][18][19].
In the parallel direction, the non-locality connects the velocity of a particle in a given point with the acceleration it has received at previous locations on the guiding center orbit, which to the lowest order in Larmor radius is parallel to the magnetic field.Unlike for the perpendicular response, the is no natural smallness parameter that can be used to expand the parallel response in the ion cyclotron frequency range, and the integral operator must be treated in full.In quasi-homogeneous plasmas, the integral operator is a convolution between the plasma response tensor and the electric wavefield.Therefore, the parallel response is often treated in Fourier space, as the convolution operator transforms into a multiplicative relation [15,16,20].In Fourier space, the plasma response tensor is a function of the wavevector, which means that each plane wave component travels with a different velocity, a phenomenon known as spatial dispersion [21].Although Fourier methods are well suited for treating the integral operator, they produce dense system matrices that are computationally expensive to invert [22], and struggle with realistic wall and antenna geometries.Furthermore, they are inefficient when describing localized small-scale phenomena such as the IBW and ICW as the discretization cannot be locally refined.To combine Fourier methods with detailed scrape-off layer simulations, it has been demonstrated that it is possible to combine two independent solutions inside and outside the plasma boundary by using mode matching techniques [23,24].
Many wave solvers neglect parallel non-local effects, with the benefit of geometrical fidelity and low computational cost [13,14,25].In practice, this can be done by using a cold plasma model (which is independent of the parallel wavenumber k ∥ ), by neglecting the poloidal field, or by guessing a value for k ∥ .In axisymmetric models, the dominating toroidal contribution to k ∥ is known exactly, whereas the poloidal contribution is not known a priori and often neglected.When the poloidal contribution is correctly accounted for, it is said to give a corresponding upshift or downshift in the parallel wavenumber.
To retain the advantages of the FEM, while including a more detailed plasma description, it is possible to add nonlocal terms to the electromagnetic wave equation iteratively.It was demonstrated in [26] that a small non-local anti-Hermitian contribution can be added iteratively to a single response tensor component in a 2D axisymmetric tokamak plasma to account for ELD in the lower hybrid range of frequencies.It has also been demonstrated that finite temperature effects, in particular Langmuir waves and ELD, can be added iteratively to a 1D finite difference model, where the initial guess is a cold plasma model [27] and the non-local contribution was computed using a kinetic description.In [27], it was also shown that without using acceleration techniques, convergence was slow and sometimes unstable.Hence, the authors applied minimum polynomial extrapolation [28] to make convergence more robust.Furthermore, it was demonstrated in [29,30] that both parallel and perpendicular non-local effects can be added iteratively to 1D axisymmetric JET-like plasmas in the ion cyclotron range of frequencies.In [30], perpendicular nonlocal effects were added to all tensor components.As such, it was possible to model mode conversion from the fast magnetosonic wave to the electrostatic IBW, as well as the propagation and damping of the IBW.Here, the non-local contribution was evaluated using wavelet transforms.To enhance convergence, the authors applied Anderson acceleration [31].
In this work, we show that it is possible to generalize this iterative method to account for non-local effects in all response tensor elements in a 2D axisymmetric tokamak plasma for ICRH applications.The non-local effects are evaluated spectrally by Fourier decomposing the electric wavefield along the magnetic flux surfaces, and the initial guess is evaluated using a local approximation of the hot plasma susceptibility tensor.The proposed iterative scheme is implemented in the full wave solver FEMIC [14] (resulting in FEMIC version 1.6) and is verified for a fast wave heating scenario in ITER.Furthermore, the iterative method will be compared to a reduced model using a local plane wave approximation (PWA).The purpose of this is twofold.First, a PWA model can be used to qualitatively verify the proposed iterative method, and second, given that the iterative method provides reliable results, it is interesting to see how well the local plane wave model performs quantitatively.The iterative model and the PWA model are compared both in terms of local power deposition and integrated power partition.
The remainder of this paper is organized as follows.In section 2, we describe the iterative method used in this work, including the evaluation of the non-local contribution as well as the local hot plasma susceptibility description.In section 3, we make a detailed study of how non-local effects change the local power deposition in a fast wave heating scenario in ITER.The iterative method is then benchmarked against the ICRF code TORIC [15] in section 4.This is followed by a discussion of the results in section 5 and conclusions in section 6.

Iterative method to account for a non-local dielectric tensor
In this section, we develop an iterative method for solving the wave equation in a 2D axisymmetric tokamak plasma with non-local response.First, in section 2.1, we describe the principles of an iterative method based on the FEM and Anderson acceleration.The plasma response is separated into a local part and a non-local dispersive contribution.Then, in section 2.2, we show how we can find a local response tensor in practice, and how to evaluate the non-local contribution using Fourier spectral decomposition.In section 2.3 we present a local plane wave model that can be used to estimate the local up-and downshift in the parallel wavenumber, and in section 2.4 we describe the local power deposition model used in this work.

Wave equation and operator splitting
The electromagnetic wave equation for a non-local plasma response is given by where J ext is the antenna current, χj is the susceptibility tensor for species j, and I is the identity tensor.Here, tilde denotes that in real space, χj is an integral operator acting on the electric field E [32].In an axisymmetric configuration, the toroidal Fourier modes decouple such that where n ϕ is the toroidal mode number and (R, ϕ, Z) are cylindrical coordinates.The toroidal mode number enters the wave equation as a parameter and the wave equation can be solved for one toroidal mode at a time.
To separate local from non-local physics, the susceptibility tensor can be split into two terms, where χ j is an local approximation of χj and δ χj is the difference between χj and χ j .Moving δ χj to the right hand side of the wave equation ( 1), we can introduce a non-local source current density to the wave equation, namely In accordance with [30], we will call this induced current the correction current density because it is correcting for errors we introduce by using a local susceptibility tensor.
As the correction current density is evaluated using the nonlocal operator δ χj , it is not natural to represent using the FEM.This can be solved by iterating between the wave equation ( 1) and the correction current (4).The electric wavefield in iteration step k would then be obtained by solving with the initial correction current density δJ (1) NL = 0.However, it was shown in [27] that a conventional fixed-point iteration scheme tends to diverge even for simple plasma models.To achieve convergence, the authors employed Anderson acceleration.Similar techniques are employed in [29,30] as well as in the present work.Adding this step to the iterative process, the electric wavefield in iteration step k is obtained by solving 2) , . .., E (1) ; where the operator A indicates Anderson acceleration and β k is a damping parameter.Note that the authors of [27] chose to iterate on the correction current density, whereas we found that iterating over the electric field provided more reliable convergence in the scenarios described in this work.Convergence can be optimized by, for example, tuning the damping parameter (which can be a function of k) and the number of stored solutions.For further details on Anderson acceleration and how it relates to other acceleration techniques, we refer to [31].

Evaluation of correction current density
The correction current density defined in (4) can for example be evaluated using poloidal Fourier decomposition, as the Fourier components are inherently non-local and are well suited to couple the current in a point r with the particles' previous acceleration in points r ′ .Using such a description, the correction current density is given by where is the Fourier representation of the susceptibility tensor χj [17][18][19], and k ⊥ and k ∥ are the components of the wavevector k that are perpendicular and parallel to the external magnetic field B, respectively.The exact implementation of χj used in this work is given below.Here, caret emphasizes dependence on k.Density and temperature profiles as well as the external magnetic field are assumed to vary slowly in space, i.e. the plasma is assumed to be quasi-homogeneous.The susceptibility tensor is evaluated in a Cartesian coordinate system (e x ,e y ,e z ), where e z = B/B, e x = k ⊥ /k ⊥ , and e y = e z × e x .The details of the coordinate transformation between this coordinate system and the cylindrical system are outlined in [33].In this coordinate system, the susceptibility tensor for a Maxwellian plasma with no net parallel flow can be represented as where Here n refers to the harmonic number, ω p,j is the plasma frequency, Ω c,j the cyclotron frequency, λ j = k 2 ⊥ ρ 2 L,j /2 the FLR parameter, ρ L,j = v ⊥,j /ω c,j the Larmor radius, v ⊥,j the perpendicular thermal velocity, v ∥,j the parallel thermal velocity, I n = I n (λ j ) the modified Bessel function of order n, Z(ζ n,j ) the plasma dispersion function evaluated at ζ n,j = (ω + nΩ c,j )/(k ∥ v ∥,j ) and ϵ j the charge sign.
For a given toroidal mode number n ϕ , the parallel wavenumber is given by where B ϕ is the toroidal magnetic field, and m θ is the poloidal mode number.In this work, we use an equal-arc-length poloidal angle θ ≡ 2π ¸dl/L, where L is the circumference of the poloidal cross section of the local flux surface.Consequently, the rotational transform ι, used to describe the local helicity of the field, can be defined as where B θ is the poloidal magnetic field.It should be noted here that due to the choice of poloidal angle, ι is not a flux function.Unlike the toroidal modes, the poloidal modes are coupled.In general, the wave exhibits a broad poloidal spectrum consisting of both positive and negative poloidal mode numbers.
The perpendicular wavenumber is approximated by the dispersion relation for the fast magnetosonic wave.To find the dispersion relation, we expand the susceptibility tensor to second order in k ⊥ and assume the parallel electric field component to be negligibly small.The dispersion equation can then be solved for k ⊥ analytically.Denoting the solution corresponding to the fast magnetosonic wave k FW = k FW (r, m θ ), we let Because the toroidal mode number enters as a parameter, the n ϕ -dependence is suppressed for brevity.The local approximation of the perpendicular wavenumber retains certain FLR effects, such as higher harmonic damping and TTMP [13,14], but does not describe mode conversion to the IBW.It is important to account for both second harmonic damping and TTMP, as they each constitute a significant portion of the total absorbed power.
To obtain a local approximation of the susceptibility tensor, we neglect the m θ contribution in (18) (and thus also the m θ dependence in k FW ).This is accurate where non-local effects are not important, or where the wave propagates perpendicular to the flux surfaces in the poloidal plane.Elsewhere, this approximation is either an under-or an overestimation of the local parallel wavenumber.This description of the susceptibility tensor is algebraic, in the sense that χ j can be evaluated algebraically in each point in space.Explicitly, the local approximation is given by This plasma description has previously been implemented in a number of wave solvers, such as LION [13] and FEMIC.
It is worth noting that if the spectrum in m θ is sufficiently wide, k ∥ can be close to zero, as can be seen from (18).If this occurs on the ion cyclotron resonance (i.e.where ω + nΩ c,j = 0), the susceptibility becomes singular, preventing the solver from converging to a physical solution.To resolve this singularity, we can expand the magnetic field along the particle orbits [34,35], using the techniques in [36] to retain the positive semi-definiteness of the imaginary part of the plasma dispersion function and thus the anti-Hermitian part of the dielectric tensor.This method introduces a toroidal broadening and is described in more detail in the appendix.
In scenarios with dominant electron damping, further issues may arise related to small values of k ∥ , this time related to the Landau response.To converge to a physical solution, it was shown in [37] that it is possible to expand the parallel wavenumber in (18) along the particle orbits.This has not been implemented in the present work, as we are only considering scenarios with dominant fast wave ion absorption.

Local plane wave model
In (21) a local approximation of the dielectric tensor was defined by setting m θ = 0.However, when the wavefield locally resembles a single plane wave, the model for the parallel wavenumber can be improved to account for the up-and downshift.It was proposed in [38] to first perform a simulation with m θ = 0 and then use derivatives of the resulting electric wavefield to estimate k.Here we propose an extension to this model that we will refer to as the plane wave approximation (PWA) model.We use derivatives of the right hand polarized wave component, E n ϕ ,− , to estimate the poloidal wavenumber, as E n ϕ ,− is not strongly effected by either the cyclotron resonance or the ion-ion hybrid layer.The estimate can be obtained from a logarithmic derivative.However, in regions with destructive interference, the wave amplitude can become small, causing the logarithmic derivative to diverge.Therefore, we here propose to estimate the wavevector in the RZ-plane as where ⟨•⟩ denotes an average over a neighborhood around r with a radius of half a wavelength, and k RZ is the projection of k onto the poloidal plane.This estimate is here used to describe the direction of k, whereas the magnitude of k is taken from the fast wave dispersion relation.Assuming that the poloidal field is small compared to the toroidal one, the wavevector in the poloidal plane is dominated by the perpendicular wavenumber, k FW .Furthermore, a weight function, w = w(r), is introduced to avoid including up-or downshift when the wave resembles an eigenmode.The parallel wavevector of the PWA model is therefore given by where we propose the following ad hoc model for the weight As an example, consider the sum of a left and a right moving wave.Let the amplitude of the right moving wave be smaller than the left moving one by factor f. It can then be shown numerically that w ≈ 1 − f for f ≲ 0.7.Thus, this model provides an average parallel wavenumber, also in the presence of wave interference.
In the context of this work, we want to use the PWA model as an alternative to the proposed iterative method.Therefore, ( 22) is applied to the right-hand polarized component of the electric field from the first iteration, E (1) n ϕ ,− .After reevaluating dielectric tensors and electric fields using k ∥ = k PWA ∥ (r), the power deposition is compared to the final iteration in the iterative method.

Power absorption
In this work, the absorbed power density by a particle species j will be evaluated as where we have used the susceptibility tensor as the power absorption kernel.It is interesting to see how the non-local contribution of the susceptibility tensor (i.e.δ χj ) affects the power absorption by particle species j, so we define the absorbed power density due to non-local absorption kernel for species j as It should be noted here that δQ j kernel is not symmetric, so we cannot expect the local power deposition to always be positive semi-definite.However, integrated absorption always is [15].
The non-local plasma response will also impact the power absorption due to changes to the electric field when iterating.Therefore, we will also define the absorbed power density due to these changes in the electric field for species j as where E (1) is the solution from the first iteration, i.e. the solution to the wave equation with δJ NL = 0.It is evident from ( 25)-( 27) that the power deposition without any nonlocal contributions is given by Q j − δQ j kernel − δQ j fields .

Effect of parallel non-locality on fast wave heating
The proposed iterative method described in section 2 has been implemented in the full wave solver FEMIC, which is based on the finite element method.The wave equation in ( 7) is discretized into triangular finite elements and solved using the RF module in COMSOL Multiphysics ® (as shown in [14]), whereas the correction current density in ( 8) is evaluated in MATLAB ® using the method outlined in section 2.2.
The correction current density is then imported to COMSOL Multiphysics ® to be used as a source term to the wave equation in the following iteration step.The impact of parallel nonlocal effects on fast wave heating is then studied in an ITERlike plasma.To analyze the effect on the different absorption modes, fundamental minority heating, ELD, TTMP, and the coupling between ELD and TTMP are all studied separately.Second harmonic tritium damping is not shown due to its qualitative similarities to the minority absorption.Although second harmonic damping is sensitive to the perpendicular wavenumber, the dispersion relation (20) does not change much for the most important mode numbers m θ .Finally, we will show the non-local effects on ELD for n ϕ = 27.For each absorption mode, we study δQ kernel and δQ fields obtained both from the proposed iterative method and the PWA model described in section 2.3 to see how well the PWA model performs compared to a more sophisticated approach.The two models are compared on a straight line starting near the bottom antenna strap.For n ϕ = 27, we compare the two models in a region without strong reflections, as this would be outside the range of validity for the PWA model.

Scenario
The effect of non-local parallel physics is studied in a fast wave heating scenario in an ITER-like plasma with a 3% 3 He minority in a DT plasma with n D = n T .The temperature and density profiles are shown in figure 1.The remainder of the model parameters are given in table 1.For n ϕ = 61 and ITER relevant temperatures, resonance layers are broad, the ion-ion hybrid layer is less pronounced compared to lower mode numbers, and there is minimal mode conversion to the ion cyclotron wave.Under these conditions, fast wave heating is straightforward to simulate using an ICRH code such as FEMIC and non-local effects can be analyzed without complications.The chosen value for n ϕ corresponds to one of the dominating toroidal mode numbers in a dipole phasing configuration (0 π 0 π).Unless stated otherwise, the results in the remainder of this section have been obtained for n ϕ = 61.The wall and antenna geometries used in this work are described in [14], and can be seen in figure 2.
Although the two poloidal antenna triplets will have a ±90 • phase shift in ITER [39], they are here simulated in phase.This simplifies the analysis and the validation of the numerical algorithm by removing the destructive interference between the top and bottom antenna triplets [14].The scenario was chosen to clearly demonstrate the qualitative effects of the  up-and downshift in the parallel wavenumber.For symmetry reasons, parallel non-local effects are not expected to have a large impact on metrics such as the total coupled power, the power partition and the power deposition profiles in this scenario [40].
The chosen scenario manifests strong single pass damping and central ion absorption, as shown in figure 2. Most of the power (65%) is absorbed by the ions, of which 90% is absorbed by the 3 He-ions and 10% is absorbed through second harmonic tritium damping.The remainder of the power is absorbed by the electrons between the antenna and the ionion hybrid layer through ELD and TTMP.
The toroidal and poloidal fields are oriented in the clockwise (negative ϕ) and counterclockwise (positive θ) directions, respectively.Consequently, for a positive n ϕ , we can expect a local upshift in the lower half of the poloidal plane, and a local downshift in the upper half.

Minority absorption for n ϕ = 61
The width of the resonance layer is approximately proportional to the Doppler shift k ∥ v ∥,j , where again v ∥,j is the thermal velocity for species j.Therefore, we expect to see a broadening of the resonance layer where the wave is upshifted (in this case, below the equatorial plane).Conversely, we expect a narrower resonance layer above the equatorial plane.This can indeed be observed in figure 3(a).In terms of local values, the contribution to the power deposition is significant, but due to the antisymmetry about the equatorial plane, the effect on integrated power deposition is small.
The parallel non-locality does not only affect where a given wavefield deposits its power in the plasma, but also the polarization of the wave.The change in polarization of the wave then impacts the power deposition in a different way than the non-locality of the power absorption kernel.A downshifted k ∥ should give a more pronounced ion-ion hybrid layer, with larger relative E + component at the hybrid layer, and smaller relative E + at the minority resonance.For an upshifted wave, we expect the opposite.This can indeed be observed when we study δQ In figure 3(c) we observe a close agreement between the FEMIC simulation and the local plane wave model.We note that the PWA model captures both the broadening of the resonance layer seen in δQ 3He kernel as well as the change in polarization represented in δQ 3He fields not only qualitatively, but also quantitatively.

Electron absorption for n ϕ = 61
An upshift in the parallel wavenumber can both enhance and reduce electron absorption, depending on the value of k ∥ v ∥,e /ω.The electron damping is more complex than the ion damping due to the different physical effects at play.The fast magnetosonic wave is damped through both ELD and TTMP.The processes are coherent, and for fast wave heating approximately twice as much power is absorbed through TTMP than ELD.ELD and TTMP are connected through the susceptibility elements χyz,e and χzy,e , and the contribution from the cross terms is negative semi-definite, approximately of the same magnitude as the TTMP term for ICRF waves in tokamaks [8].All three contributions to the electron absorption have different dependencies on k ∥ v ∥,e /ω [17][18][19].Therefore, we will study each absorption mode separately.
For the electron Landau damping, we see in figure 4(a) that an upshifted parallel wavenumber leads to a negative contribution from the non-local absorption kernel, δQ ELD kernel .Conversely, a downshifted parallel wavenumber yields a positive contribution.As for the effect of the iterated electric fields, δQ ELD fields , we see the opposite behavior (see figure 4(d)), although the contribution is not large enough to fully cancel the effect of the non-local absorption kernel.The non-local contribution from the cross terms is structurally very similar to the contribution from the ELD, although with opposite signs (see figures 4(c) and (f )).This is explained by the similar (but not identical) functional dependence on k ∥ v ∥,e /ω.
Regarding the power deposition due to TTMP, we see in figure 4(b) that δQ TTMP kernel vanishes at R = 7.5 m.This is due to passing a maximum in the zeroth harmonic term of the susceptibility tensor component χyy,e .Again, we note an antisymmetry about the equatorial plane.We also observe in figure 4(e) that the non-local effect from the iterated electric fields changes sign at the ion-ion hybrid layer.This is due to a similar polarization change that was observed for the ion absorption in the previous section, but now for the E ycomponent.
The PWA model produces a fairly close agreement with the iterative model.This is shown for the different terms in figures 4(g)-(i).The largest differences are observed between R = 6.75 m and R = 7.25 m due to the reflection of the fast magnetosonic wave off the ion-ion hybrid layer, resulting in an interference pattern.

Electron Landau damping for n ϕ = 27
Given the agreement between the plane wave model and the implemented iterative method, we can use (23) to evaluate the relative shift in k ∥ for n ϕ = 27 and n ϕ = 61, respectively.The relative upshift varies over the poloidal plane depending on how well the wave fronts are aligned with the flux surfaces.We find that for n ϕ = 27, the relative upand downshift is approximately twice as large as for n ϕ = 61.This can be confirmed using (18) for a given poloidal m θ spectrum.The larger shift can also be seen in the magnitudes of δQ ELD kernel and δQ ELD fields in figures 5(a) and (b).The non-local contributions are more than twice as large compared to the corresponding terms for n ϕ = 61 (figures 4(a) and (d)).
In figure 5(b) we also see that the reflection of the fast wave off the ion-ion hybrid layer produces a wave pattern in δQ ELD fields near the equatorial plane.The incoming wave from the antenna is upshifted, whereas the reflected wave is downshifted.This alters the corresponding interference pattern in the electric fields and electron absorption.Interference is not well modelled by a plane wave.We therefore make the comparison with the PWA model below the region of strong interference.Despite the large upshift, we see in figure 5(c) that the PWA model still produces good agreement with the iterative method, although it seems that the PWA model predicts some interference where the iterative method predicts none.Similar comparisons can be made for the other absorption modes as well.
It is noteworthy that whereas the up-and downshift is more perturbative in nature for n ϕ = 61, the effect for n ϕ = 27 is of zeroth order locally, completely changing the local power deposition.Ion absorption for n ϕ = 27 is not shown due to its qualitative similarities to n ϕ = 61.The main difference is a smaller Doppler broadening of the resonance for n ϕ = 27 due to smaller k ∥ .

Convergence
The relative residual norm for iteration step k is defined as L 2 has been rescaled by its L 2 -norm to ensure that the left-handed and in particular the parallel components are not neglected when the residual is evaluated.The iteration procedure is deemed to have converged when the relative residual norm is smaller than a set value, in this case 10 −3 .The proposed numerical algorithm shows good convergence properties, as shown in figure 6.The rate of convergence is dependent on the toroidal mode number n ϕ , as a lower toroidal mode number such as n ϕ = 27 leads to a larger relative up-or downshift in the parallel wavenumber, requiring more iterations.It is also possible that increased mode conversion to the ICW for lower mode numbers slows down convergence.However, the most likely source to the slow convergence is small amounts of unresolved mode conversion near the magnetic axis, which can be seen in figures 5(a) and (b).This introduces numerical noise to the iteration process, making it more difficult to converge.
To reduce numerical noise, is possible to filter the correction current density after each iteration.In this work, we have used total variation denoising (TVD) to filter out numerical noise from the finite element solution.It is in some cases possible to reduce the number of iterations significantly, as shown in figure 6, with little impact on the resulting wavefield and power deposition.For n ϕ = 27 we observe an improvement of almost a factor of two, whereas for n ϕ = 61 the difference is marginal.There is further room for optimization of the convergence, and possible parameters to explore are the TVD filtering parameter, the damping parameter, and the number of stored solutions in the Anderson acceleration.Such optimization is however outside the scope of this paper and may be addressed in a future publication.Convergence is much slower for lower values of n ϕ , but can be sped up by filtering the correction current density using TVD filtering (dashed curves).

Comparison between FEMIC and TORIC
To quantitatively validate the iterative method proposed in this work, FEMIC is benchmarked against the full wave ICRF code TORIC.TORIC uses a finite element discretization in the radial direction and Fourier methods in the poloidal and toroidal directions.The spectral treatment in the poloidal direction includes non-local effects natively, such as the up-and downshift in the parallel wavenumber.TORIC uses a second order FLR expansion in the perpendicular direction to account for higher harmonic damping and mode conversion to the IBW [41].Unlike in FEMIC, perpendicular non-local effects (FLR effects) are treated with differential operators.FEMIC's local treatment describes higher harmonic damping well but does not model mode conversion to the IBW.
In this exercise, we compare three different models for parallel response in FEMIC to TORIC.First, we study FEMIC without any iterative additions to estimate how important nonlocal effects are in this particular scenario.We then apply the proposed iterative method, and finally, we also consider the PWA introduced in section 2.3.
The models are benchmarked by comparing their predicted power partition.Here we use the same scenario as in section 3, but sweep over minority concentrations from 0% to 4% in halfinteger steps.To avoid any of the antisymmetries observed in section 3, and thus enhance the non-local effects on the power partition, the whole antenna is located below the equatorial plane.This antenna model is similar to using only the bottom triplet of the ITER antenna.To represent the same geometry in the two codes and thus facilitate the comparison between them, a simplified wall and antenna geometry has been implemented in FEMIC.Furthermore, we set the toroidal mode number to n ϕ = −27 to get a downshift in the parallel wavenumber.This comparison is similar to the exercises done in [14,42].

Power partition
In the pure second harmonic heating scenario, TORIC predicts approximately 5% of the power to be mode converted into the IBW.The IBW is mostly absorbed by the electrons on the high field side.When a small minority concentration is introduced, the power converted to the IBW rapidly shrinks to less than 1% of the total deposited power.As FEMIC does not support IBW physics, the comparison of the power partition predicted by the two codes is facilitated by neglecting the IBW electron absorption in TORIC.To account for the removed IBW absorption, the power partition in TORIC is renormalized.The result is shown in figure 7.
When non-local effects are not accounted for, FEMIC systematically overpredicts the electron absorption compared to TORIC for all minority concentrations, ranging from 15% in the pure second harmonic tritium heating scenario to 5% at 4% 3 He concentration.This is not surprising considering the large downshift in the parallel wavenumber below the equatorial plane.Consequently, FEMIC also underpredicts the fraction of power absorbed by the ions.
Including parallel non-local effects iteratively, we note a close agreement between FEMIC and TORIC for all particle species, in particular for lower minority concentrations.The difference is the largest for 4% 3 He concentration but is even then only a few percent.It is not surprising that we observe slightly poorer agreement for higher concentrations.FLReffects tend to be more important here, and TORIC has a more advanced description of the perpendicular non-local physics.
The local plane wave approximation of k ∥ shown in blue reproduces the good agreement of the iterative method for all minority concentrations, including concentrations below 1%, where single pass damping is weak.Here, the wave has no clear beam shape and is poorly approximated by a plane wave and we can expect the width of the poloidal spectrum to be important for wave polarization and absorption.Still, the difference in electron absorption between FEMIC and TORIC is small.
All four models predict maximum absorption by the 3 Heions near 2.5% 3 He concentration.As the minority concentration increases, the ion-ion hybrid resonance moves away from the cyclotron resonance.Consequently, the overlap between the cyclotron resonance layer and ion-ion hybrid layer decreases, and the favorable polarization is not exploited.

Discussion
The main advantage of the implemented iterative method is the possibility to model the complex wall and antenna geometries in high fidelity thanks to the memory efficient finite element method, while still able to account for nonlocal effects, which are added iteratively.physics is only accounted for in the plasma core and is not expected to be significant in the SOL due to the low plasma temperatures there, saving computational resources.Another advantage of using a FEM solver is the possibility for mesh optimization, such as local mesh refinement for short wavelength wave modes.Using the FEM in the whole poloidal plane, combined with Fourier decomposition to account for parallel nonlocal effects, the proposed iterative method is less memory intensive than many other comparable methods at the expense of longer computational time.In principle, this is a desirable trade-off making it possible to run the code on a workstation, rather than a cluster.This applies also when modelling large machines, such as ITER and the EU-DEMO.Furthermore, fixed point iteration methods (and sometimes Anderson acceleration in particular) are being employed in other codes to couple the wave solver to a Fokker-Planck solver [43], or to apply rectified sheath boundary conditions [44,45].FEM solvers work well in this type of framework.Each wave solution typically requires a few minutes of computational time, and in principle, non-local effects, sheath physics and non-linear Fokker-Planck physics can all be included in each iterative step, reducing the total computation time compared to a more expensive wave solver.Unless the separate effects are strongly coupled, the convergence and total computation time should be mostly determined by the physical effect that converges the slowest.
The reduced local plane wave model, which was used to validate the proposed iterative method, produces surprisingly good agreement not only in terms of power partition as shown in the benchmark in section 4, but also in local power density as shown in section 3.This good agreement is not expected to be maintained in scenarios with poor single pass damping, especially in smaller machines.When single pass damping is low, the reflected and refracted waves develop a complex poloidal mode spectrum, and a local plane wave model is then not sufficient to represent the up-and downshift in the parallel wavenumber.
An attempt was made to use the PWA model to improve the convergence rate by providing an improved starting guess, but the error norm was only reduced in the first few iterations.The local plane wave model prescribes a k RZ mainly in the negative R-direction.Therefore, the reflected portion of the wave gets the opposite shift in the parallel wavenumber compared to what it should have.The iterative method then has to compensate for such errors, negating any gains from the improved initial guess.In fact, in some scenarios, the rate of convergence was reduced when using the PWA model as an initial guess.Despite this, the PWA model remains a useful first order correction for parallel non-local effects, especially in large machines.The reason for this is twofold.First, the model is more accurate in large machines with strong single-pass damping, such ITER and the EU-DEMO.Second, methods that take non-local effects into account tend to be computationally expensive, increasing the potential gain by using memory efficient methods such as the FEM.In particular, 3D FEM models can benefit from such reduced models [25].
In terms of relative importance of the various contributions the total deposited power, it is interesting to note that the quantities δQ j kernel and δQ j fields are of similar magnitude for a given species j.This means that accounting for nonlocal effects in the absorption kernel and in the wave propagation is equally important.Furthermore, comparing δQ j kernel and δQ j fields (figures 3 and 4) to Q j (figure 2), we note that the relative change in absorption due to non-local effects is similar for all absorption modes.This implies that non-locality in all tensor components is equally important, perhaps with the exception of the xz-and zx-components due to their small magnitude.All other components significantly contribute to the power deposition.
A final point should be made on the importance of the upand downshift on the parallel wavenumber.For dipole phasing, the non-local effects are perturbative in nature, slightly changing the resonance width for ions and the absorbed power density for the electrons.For toroidal phasings that have a larger portion of the spectrum in the lower range of toroidal mode number, such as [0 π π 0], [0 0 π π] or current drive phasing, we can expect a larger effect.The nonlocal contribution to the local power density may then be of zeroth order, significantly changing the total absorbed power density.
The comparison with TORIC produced good agreement in general, but in particular for lower minority concentrations where FLR effects are less important.At 3%-4% 3 He concentration it may be important to include FLR physics.This may be difficult to implement in FEMIC.Many of the terms of interest are divergence-like terms, which are difficult to evaluate in the RF module in COMSOL Multiphysics® as the curl-conforming vector elements used to discretize the electric field are not well suited for this type of operations.Here other basis functions or other formulations may need to be considered.It is also possible that the discrepancy is due to different representations of the parallel physics.

Conclusions
We have implemented an iterative method to include parallel non-local physics in the FEM based ICRF code FEMIC.The non-local contribution was evaluated using Fourier decomposition in the poloidal angle coordinate.The new version of FEMIC was applied to a fast wave minority heating scenario in ITER and has been benchmarked against the ICRF code TORIC with good agreement.To further verify and explain the change in wave polarization and local power absorption due to the inclusion of parallel non-local effects, a reduced local plane model was developed and compared to the simulation results from FEMIC.We showed that we can accurately represent the broadening of the resonance layer, the change in polarization due to the non-local effects on the ion-ion hybrid layer, as well as capture the impact of non-local effects on the different electron absorption modes.We also show that using a non-local absorption kernel is equally important as using a non-local dielectric tensor, for both ion and electron absorption.
One of the main advantages of the proposed numerical scheme is that the finite element method in 2D axisymmetric configurations (or 3D) is well suited to be applied to complex geometries.Wave coupling and propagation in the scrape off layer and at the plasma edge are well represented by the finite element method without including non-local effects, as the temperature is low enough that non-local effects become negligible.Non-local contributions only need to be added in the plasma core, where the temperature is high.Furthermore there is no need to invert a large dense system matrix, as the non-local effects are added iteratively.In principle we can add more physical effects in the iteration procedure, such as nonlinear boundary conditions, or coupling the wave solver to a Fokker-Planck solver.
The current implementation lacks certain finite Larmor radius effects, such as mode conversion to the ion Bernstein wave (IBW).This could be included by expanding the response tensor in k ⊥ and giving higher order terms a finite element representation so that they can be evaluated in a FEM solver such as COMSOL Multiphysics®, although care needs to be taken to use appropriate basis functions.The mode conversion to the IBW is sensitive to the up-and downshift in the parallel wavenumber [15], effects that in principle can be included using the method in the present work.

Figure 1 .
Figure 1.Profiles of electron density (black), electron temperature (dashed orange) and ion temperature (dash-dotted blue) as a function of major radius in the equatorial plane.

Figure 2 .
Figure 2. Real part of the left hand circularly polarized component of electric field, (a) E+, and absorbed power density by (b) 3 He-ions and (c) electrons.Minority and ion-ion hybrid resonances are indicated by dashed and dotted vertical lines, respectively.

Figure 3 .
Figure 3. Contribution to the absorbed power density by 3 He-ions due to non-local physics in (a) the absorption kernel (δQ 3 He kernel ), and (b) the susceptibility tensor (δQ 3 He fields ) as predicted by FEMIC in units of mWm −3 /W.Minority and ion-ion hybrid resonances are indicated by dashed and dotted vertical lines, respectively.In subfigure (c), the iterative method (blue and black) is compared to the local plane wave model (red and green) along the oblique line (black, figures (a) and (b)).Here n ϕ was set to 61.

Figure 4 .
Figure 4. Contribution to the absorbed power density through ELD (first column), TTMP (second column) and cross terms (third column) due to non-local physics in the absorption kernel (δQ kernel ) (a)-(c), and susceptibility tensor (δQ fields ) (d)-(f ) as predicted by FEMIC in units of mWm −3 /W.Minority and ion-ion hybrid resonances are indicated by dashed and dotted vertical lines, respectively.In subfigures (g)-(i), the iterative method (blue and black) is compared to the local plane wave model (red and green) along the oblique line (black, figures (a)-(f )).Here n ϕ was set to 61.

Figure 5 .
Figure 5. Contribution to the absorbed power density through ELD due to non-local physics in (a) the absorption kernel (δQ ELD kernel ), and (b) the susceptibility tensor (δQ ELD fields ) as predicted by FEMIC in units of mWm −3 W −1 .Minority and ion-ion hybrid resonances are indicated by dashed and dotted vertical lines, respectively.In subfigure (c), the iterative method (blue and black) is compared to the local plane wave model (red and green) along the oblique line (black, figures (a) and (b)).Here n ϕ was set to 27.

Figure 6 .
Figure 6.Relative residual norm as a function of iteration number.Convergence is much slower for lower values of n ϕ , but can be sped up by filtering the correction current density using TVD filtering (dashed curves).

Table 1 .
Summary of model parameters.