Iterative addition of finite Larmor radius effects to finite element models using wavelet decomposition

Modeling the propagation and damping of electromagnetic waves in a hot magnetized plasma is difficult due to spatial dispersion. In such media, the dielectric response becomes non-local and the wave equation an integro-differential equation. In the application of RF heating and current drive in tokamak plasmas, the finite Larmor radius (FLR) causes spatial dispersion, which gives rise to physical phenomena such as higher harmonic ion cyclotron damping and mode conversion to electrostatic waves. In this paper, a new numerical method based on an iterative wavelet finite element scheme is presented, which is suitable for adding non-local effects to the wave equation by iterations. To verify the method, we apply it to a case of one-dimensional fast wave heating at the second harmonic ion cyclotron resonance, and study mode conversion to ion Bernstein waves (IBW) in a toroidal plasma. Comparison with a local (truncated FLR) model showed good agreement in general. The observed difference is in the damping of the IBW, where the proposed method predicts stronger damping on the IBW.


Introduction
Modeling the propagation and damping of electromagnetic waves in a hot magnetized plasma is difficult due to spatial dispersion [1][2][3][4]. Spatial dispersion occurs when natural length scales are present in a medium. As a consequence, the dielectric response becomes non-local and the wave equation an integro-differential equation that is difficult to solve with common numerical methods [5][6][7].
In the application of RF heating and current drive in tokamaks [8], accounting for spatial dispersion is important in modeling codes. Parallel spatial dispersion is caused by the parallel velocity of the particles, and is important for Landau damping and the Doppler broadening of the ion cyclotron resonances. In the presence of a poloidal magnetic field, parallel spatial dispersion makes the dielectric tensor an integral operator [7]. In the perpendicular direction, the finite Larmor radius (FLR) of the gyrating particles causes spatial dispersion. FLR effects give rise to higher harmonic ion cyclotron damping, transit-time magnetic pumping (TTMP) and mode conversion [1,2,9,10].
Common numerical tools for solving the wave equation include the finite element (FE) method and Fourier spectral methods. The FE method is efficient in solving the wave equation and can be applied on complex geometries. The main disadvantage of the FE method is that it cannot account for the non-local character of the dielectric response. Fourier spectral methods can describe non-local effects by using Fourier bases, however, these methods are time consuming and require a large amount of computer memory.
Much effort has been devoted to develop modeling codes that can simulate the propagation and damping of RF waves in hot magnetized plasmas [11][12][13][14][15][16][17][18]. A common approach is to solve the wave equation in two dimensions and take advantage of axisymmetry in the toroidal direction. The advantage of axisymmetry is that the electric field can be Fourier decomposed toroidally and simulations can be performed for a single toroidal mode number at a time. In the radial and poloidal directions (or poloidal plane), the choice of discretization method varies between different modeling codes. The simplest approach is to apply a FE discretization in the poloidal plane. This method relies on algebraic approximations that remove the integral operators in the dielectric tensor, so that the wave equation becomes a partial differential equation (e.g. by neglecting the poloidal contribution to the parallel wave number) [16,18]. One approach to include nonlocal effects due to parallel dispersion is to Fourier decompose in the poloidal direction and use a FE discretization in the radial direction [11,12]. For this type of numerical scheme, FLR effects are typically included by expanding the dielectric tensor in the Larmor radius, where the first few terms are kept. Such approximations are sufficient to describe e.g. higher harmonic damping and mode conversion from the fast magnetosonic wave (or fast wave) to electrostatic waves, but are only valid when the Larmor radius of the gyrating ions is relatively small. Fourier spectral methods can solve the integro-differential wave equation without imposing any approximations [17]. However, these methods require massive computing power and are difficult to apply outside the plasma domain.
Research in development of new numerical schemes to include non-local effects using the FE method has been made and is available in the literature. In [19,20], Sauter et al developed a method for solving the one-dimensional wave equation using a FE scheme, where FLR effects are included by employing a gyrokinetic approach. In [21], a numerical method called Gabor element method was developed, which was used to study a case of mode conversion from fast to slow waves. In this numerical method, the electric field is expanded into Gabor functions (note that Gabor functions are similar to Morlet wavelets used in the present work), which results in sparse matrices. In [22], Meneghini et al developed a twodimensional FE model for lower hybrid waves in tokamaks. This method is based on a cold plasma model, where hot plasma effects are added by iterations to the most important dielectric tensor elements. A similar iterative algorithm in one dimension was developed by Green et al [23], where parallel non-local effects are evaluated by calculating particle orbits backward in time and added by iterations to the wave equation. Recently, an iterative wavelet FE scheme has been developed that is suitable for modeling RF heating of tokamak plasmas [24]. This numerical scheme was tested for parallel dispersion in one dimension, where non-local effects were evaluated using wavelets and added to the wave equation by iterations.
In this paper, the main objective is to show that the iterative wavelet FE scheme can account for spatial dispersion caused by FLR effects in one dimension. We split the susceptibility operator into an algebraic susceptibility tensor and an induced current density term (were the latter quantity contains the integral operators). The wave equation is solved with a FE scheme using the algebraic susceptibility tensor, where the induced current density is represented as an inhomogeneous term. The contribution of the induced current density is added by using a fixed point iteration scheme based on Anderson Acceleration [25]. The induced current density is assumed to be zero for the first iteration. For subsequent iterations, the induced current density is evaluated using solutions from previous iterations. To account for the nonlocal character, the induced current density is evaluated in wavelet space [24,[26][27][28][29]. In this work, the Morlet wavelet was used to perform the wavelet decomposition. Morlet wavelets have several attractive features, such as having good convergence properties in the wavelet transform as well as being localized in both configuration and frequency space [30][31][32].
As an example, we will apply this method to a case of fast wave heating at the second harmonic ion cyclotron resonance, and study mode conversion to ion Bernstein waves (IBW) for different toroidal mode numbers. In this model, we use the hot susceptibility tensor for Maxwellian plasmas, including FLR effects. We show that the solution converges and that it is consistent with the hot plasma dispersion relation. The solution is also compared with the results from the TOMCAT code [33], which is a one-dimensional ion cyclotron resonance heating code.
The paper is structured as follows: the wave equation and susceptibility tensor are described in section 2; the continuous wavelet transform and the wavelet representation of the susceptibility tensor are described in section 3; the simulation results in section 4, and the discussion and conclusions in sections 5 and 6, respectively.

Wave equation and susceptibility tensor
In this section, we develop a one-dimensional model for ion cyclotron resonance heating in a toroidal plasma. We derive an expression for the wave equation that is suitable to be solved with the iterative wavelet FE scheme. Furthermore, we describe the boundary conditions and give expressions for the susceptibility model. The last sections describe the hot plasma dispersion relation and the energy equation.

Wave equation
The wave equation is given in a cylindrical coordinate system ( ) f e e e , , R Z , where e R is the major radius direction, e f the toroidal direction and e Z the vertical direction. We assume that the plasma response and geometry are axisymmetric (i.e. J E . 5 ind j j j 0 To solve equation (4), we iterate over δJ ind in a fixed point iteration scheme using Anderson Acceleration [24,25]. Let k be the iteration index and E (0) =0 the initial solution. A solution to the wave equation can then be found by solving the following equations, Thus, equation (7) can be solved using the FE method. To evaluate equation (6) and account for the non-local character of the equation, we use a method based on Morlet wavelet decomposition, which is described in section 3.

Boundary conditions to the wave equation
In this model, we assume that the domain is confined by perfectly conducting walls. Thus, the boundary condition applied is where n is the unit normal vector to the wall. The antenna is located inside the domain and is represented as an internal boundary, where we assign a prescribed surface current density. In cylindrical coordinates, we let the antenna location be described by R ant . The antenna surface current density J ant is then given by where J 0 is a complex constant with the unit A m −1

Hot susceptibility tensor
In this work, we use a quasi-homogeneous approximation of the susceptibility tensor, where we assume a slow variation in the density, temperature and magnetic field. The susceptibility tensor is given in a Cartesian coordinate system ( ) e e e , , x y z , where the magnetic field is aligned in the positive e z direction and the wave vector is in the xz plane (i.e. k y =0). Denoting Fourier transformed variables using the hat symbol, the susceptibility tensor for a hot Maxwellian plasma is given by [ where the elements are given bŷ the FLR parameter, ρ L,j is the Larmor radius, v j the thermal velocity, I n =I n (λ j ) the modified Bessel functions of order n, and Z(ζ n,j ) the plasma dispersion function evaluated at z = n j , . The prime denotes derivative. The magnetic field in tokamaks can be decomposed into toroidal and poloidal components. In this work, we neglect the contribution of the poloidal magnetic field in order to focus on FLR effects only. Such an assumption is reasonable, given that the poloidal magnetic field in general is much weaker than the toroidal magnetic field. This assumption also simplifies the transformation matrix from the Cartesian coordinate system to the cylindrical coordinate system, and removes the integral operators due to parallel dispersion (the proposed numerical scheme was tested for parallel dispersion in [24]).
Since the wave equation is described in a cylindrical coordinate system, the susceptibility tensor must be transformed. Let e x be parallel to e R , e f antiparallel to e z , and e y parallel to e Z . Under this assumptions, the transformation matrix becomes ⎡ and the transformed susceptibility tensor is given bŷ Without a poloidal magnetic field, the parallel wave number is antiparallel to the toroidal wave number. Hence, the parallel wave number becomes

Algebraic susceptibility tensor
The algebraic susceptibility tensor used here is based on the hot susceptibility tensor from previous subsection, but evaluated for the fast wave. This requires that the perpendicular wave number of the fast wave is known in advance [18]. The perpendicular wave number of the fast wave can be obtained from the following approximative dispersion relation ⎡ , . An approximative algebraic susceptibility model can then be obtained by evaluating the hot susceptibility tensor for the fast wave, i.e.
The main advantage of using this susceptibility tensor instead of a cold plasma model is that higher harmonic cyclotron damping, TTMP and Landau damping are described quite well for the fast wave.

Hot plasma dispersion equation
The local dispersion equation is obtained by setting the determinant of the Fourier transformed wave equation to zero (equation (1)). Using the Cartesian coordinate system, the dispersion equation is given by whereĉ = + å K I j j , n x =ck x /ω and n z =ck z /ω. The roots were calculated as the intersection of the contours Re(D)=0 and Im(D)=0 in the complex k x plane.

Energy equation
For time harmonic waves in spatially dispersive media, the energy equation has the form where S is the Poynting flux, T the kinetic flux and P abs is the local absorption. The Poynting flux is given by Expressions for the kinetic flux and the local absorption in inhomogeneous spatially dispersive media are given in [6,34,35]. In these expressions, the electric field is Fourier expanded and the kinetic flux and local absorption are calculated by using the absorption kernel. Due to the complexity of the theory, an equivalent expression based on Morlet wavelets has not been derived yet. In this work, we will use approximative expressions to evaluate the kinetic flux and the local absorption, which can be written as H andĉ j A are the Hermitian and anti-Hermitian parts of the susceptibility tensor, respectively. In these expressions, are evaluated using wavelet decomposition (i.e. non-locally). The electric field component E * on the left side of these factors is evaluated locally. Therefore, the evaluation becomes non-symmetric, and the results will only be approximative. However, for the fast wave, T≈0, so that ·  » -P P abs , and as such, this approximation should be reasonable. For the IBW, » P 0, so that ·  » -P T abs . Thus, the divergence of the kinetix flux of the IBW gives the local absorption.

Wavelets and the susceptibility tensor
In this section, we will start with a short description of the Morlet wavelet and the continuous wavelet transform, followed by a derivation of the susceptibility tensor in wavelet space. We also describe the wavelet grid, matched layers and how to discretize the wavelet scale parameter.

Morlet wavelet
Let ψ(x) be the mother wavelet, which is a complex and continuously differentiable function. To generate a set of wavelets (or daughter wavelets) from the mother wavelet, we introduce the wavelet scale parameter a and the translation parameter b. The set of wavelets is then given by [30][31][32] In this paper, we use the normalized complex Morlet wavelet that is defined as where σ=6 is a reference wave number and The central wave number of each Morlet wavelet is given by k=σ/a. The Fourier transformed Morlet wavelet iŝ

Continuous wavelet transform
The continuous wavelet transform of a continuously differ- The inverse continuous wavelet transform returns the original function, given by where the admissibility constant is

Wavelet representation of the susceptibility tensor
As described in section 2.1, the considered plasma is assumed to be inhomogeneous in the e x direction only, and the wave vector to lie in the xz plane, i.e. k y =0. Under these circumstances, we can Fourier transform the susceptibility operator acting on the electric field in the e x direction as follows [6]ˆ( Since the plasma is homogeneous and periodic in the toroidal direction, the susceptibility can be evaluated for one toroidal mode number at a time. We hereby drop the k z dependence for simplicity in writing. By rewriting the electric field as a Fourier transform, we obtainˆ( We can now wavelet decompose the electric field using equation (31). Furthermore, we note that the integral in x′ is the Fourier transform of the Morlet wavelet. The result iŝ Next, we Taylor expand the susceptibility at k x =σ/a here p is the expansion order. Performing the integral over k x requires some technical steps presented in appendix. The result is where c j ab , is the susceptibility represented in wavelet space given by where He p are the probabilists Hermite polynomials [38].
Substituting into equation (33) gives the following result, Equation (39) implies that the non-local susceptibility response can be evaluated by multiplying the susceptibility and electric field in wavelet space, and then performing the inverse wavelet transform.

Evaluation of the correction current
To evaluate the correction current, we use equation (39) in (6), which gives is the wavelet transform of the electric field.

The wavelet grid and matched layers
It is not necessary to evaluate δJ ind over the whole domain used in the FE discretization. The FLR effects are important where the plasma is hot, whereas in regions with low temperature the FLR effects are expected to be small. It is therefore convenient to define a smaller grid for the evaluation of δJ ind in order to reduce the computational cost. This grid is called the wavelet grid.
Another complication occurs when a finite wavelet grid is used in the continuous wavelet transform, which is defined on an infinite interval. The introduction of a finite grid results in edge artifacts when performing the wavelet transform. These artifacts occur when the wavelets are near the boundaries, so that information outside the finite interval is required. In many applications, this is solved to some degree by using extension methods, which include zero padding, symmetric and periodic extensions [39,40]. In this work, we solve this issue by using a technique called matched layers at the boundaries of the wavelet grid [24,26]. This is basically an extension method where the function of interest is extended by matching it with a Morlet wavelet with similar wave number (see figure 1). The method ensures that the function becomes continuous over the boundary and that all edge artifacts in the wavelet transform fall within the matched layers. The evaluation of δJ ind is therefore performed in the extended wavelet grid. The matched layers are then removed on the correction current before solving the wave equation.

Scale parameter and admissibility constant
The scale parameter a is assumed to be continuous in the wavelet transform. In order to perform the calculations numerically, the scale parameter must be discretized. The calculations in this paper are based on a discretized scale parameter given by where m is an integer. The constant C ψ in equation (31) depends on the discretization and have been numerically calculated to be C ψ ≈1.1598.

Fast wave heating and mode conversion
To test the proposed numerical scheme, the code FEMIC1D was developed in MATLAB ® . In FEMIC1D, the wave equation was solved using the RF module in COMSOL Multiphysics ® , while the evaluation of equation (6) was performed in MATLAB ® . The interface between COMSOL Multiphysics ® and MATLAB ® was handled using the LiveLink TM for MATLAB ® .
The model parameters, summarized in table 1, were chosen to mimick a JET-like scenario. We considered a pure hydrogen plasma, where the antenna frequency and magnetic field were tuned for second harmonic heating of hydrogen ions near the magnetic axis. We included up to 8 harmonics in the susceptibility tensor (equation (10)), and four terms in the wavelet representation of the susceptibility tensor (equation (38)). The toroidal magnetic field was given by where B 0 is the on-axis magnetic field and R 0 is the major radius. The density and temperature profiles are shown in figure 2. For simplicity, flat profiles were assumed in the plasma with a smooth transition over the pedestal to the lowdensity domain (or scrape-off-layer). The pedestal was described by an mtanh model [18]. The electron density in the low-density domain was n SOL =4. . To solve the wave equation using the FE method, a grid size of 0.5 mm was applied in the plasma domain and 5.0 mm in the low-density domains. Second order elements were used. The wavelet grid was defined between   R 1.96 3.84 m, with a grid step of 0.63 mm. The scale parameter was chosen so that wave numbers ranging between 20 and 1536 m −1 were resolved.
The model with n f =15 required a wavelet grid with finer resolution. The step was set to 0.54 mm and modes with wave numbers between 20 and 1827 m −1 were resolved.

Dispersion relation
The dispersion relation for n f =20 is shown in figure 3(a). The fast wave branch has wave numbers ranging between   k 30 6 5 x in the plasma, and is evanescent in the lowdensity regions. Its imaginary component indicates that the second harmonic ion cyclotron resonance is located at R=3.1 m. A mode conversion zone can be identified at the resonance, where the dispersion relations of the fast wave and the IBW intersect. The IBW propagates in the negative major radius direction and becomes strongly damped near R≈ 2.2 m. Note that the imaginary component of k x for the IBW is negative, indicating that the wave is backward propagating  A third branch is observed near R≈3.8 m, which is an IBW associated with the third harmonic ion cyclotron resonance. However, this wave was not excited in the plasma and will not be discussed further. Figure 3(b) shows the dispersion relations for the IBW for the toroidal mode numbers 15, 20 and 27. The real parts of the dispersion relations are only weakly dependent on the toroidal mode number. In contrast, the imaginary part defining the damping location is sensitive to the toroidal mode number. The lower the toroidal mode number, the longer distance the IBW can propagate from the second harmonic cyclotron resonance before being strongly damped.

Electric field for n f =20
The solution to the wave equation for n f =20 as predicted by FEMIC1D is shown in figure 4. The incoming fast wave couples to the plasma and propagates in the negative major radius direction. The wave becomes damped when crossing the ion cyclotron resonance near R=3.1 m, which we identify as a reduction of the wave amplitude (there is also a weak change in the amplitude due to the change in the toroidal magnetic field). The fast wave is reflected at the pedestal near R≈2.0 m and crosses the second harmonic cyclotron resonance a second time.
The fast wave is partly converted to IBWs when crossing the second harmonic cyclotron resonance. The presence of the IBW is only observed in the E R and E f components as a short wavelength mode, which is a typical feature of electrostatic waves [41]. Figure 5 shows the propagation and damping of the IBW, which was obtained by filtering the electric field using wavelets. The IBW is born near the second harmonic cyclotron resonance and propagates in the negative major radius direction. We note that the E f component of the IBW increases when propagating away from the resonance. On the other hand, the E R component increases rapidly at the resonance and varies weakly when propagating away from the resonance. Thus, since electron Landau damping is quadratic in the amplitude of E f , the damping of the IBW becomes stronger. When the wave crosses R≈2.4 m, it becomes strongly damped and vanishes. This behavior is consistent with the hot plasma dispersion relation in figure 3.
The amplitude of the wavelet coefficients of the E R component is shown in figure 6. The incoming fast wave branch is seen in the lower half of the spectrogram, which corresponds to waves with negative wave number. The ion cyclotron damping at R≈3.1 m can be identified as a reduction of the amplitude of the wavelet coefficients. At the plasma edge near R≈2.0 m, we can observe that the fast wave is reflected (due to the cutoff) and reappears in the upper half of the spectrogram. The IBW branch is here only seen in the upper half, which corresponds to waves with positive wave numbers. However, since the IBW propagates in the negative major radius direction (i.e. negative group velocity), this is yet another indication that the wave is backward propagating.
In figure 7, the electric field components as predicted by TOMCAT are shown. By comparing the E Z component in figure 7(b) with figure 4(c), we find good qualitative agreement between FEMIC1D and TOMCAT, which indicates that the description of the fast wave is similar. By comparing the E R components in figures 7(a) and 4(a), we observe good agreement for the fast wave only. The main difference is the propagation and damping of the IBW. The characteristic of the IBW as predicted by TOMCAT is shown in figure 8. Near the second harmonic resonance, we observe a similar behavior between the two codes. When propagating away from the second harmonic resonance, the amplitude in E R stays constant down to the plasma edge. As for the E f component, we note that its amplitude grows as it propagates away from the resonance. However, despite that E f is growing, the electron Landau damping does not damp out the wave. We also note that the amplitude of the IBW near the resonance in TOM-CAT is about a factor two larger than in FEMIC1D, indicating a larger fraction of mode conversion. A possible explanation is that the interference patterns caused by the incoming and reflected fast waves are different in FEMIC1D and TOMCAT. As was shown in [42], the mode conversion rate is sensitive to these interference patterns.

Power absorption
The total power absorption P abs and absorption profiles were calculated assuming unit current density in the antenna strap. Since we have a one-dimensional axisymmetric model where the solution is invariant in the e Z direction, the total power absorption was estimated by integrating the power absorption density in the Rf plane (i.e. P abs has the unit W m −1 ). The results are summarized in table 2 and show that the absorbed power increases with decreasing toroidal mode number. This behavior is consistent with both theoretical predictions and experiments [18,43,44]. Furthermore, both FEMIC1D and TOMCAT predict that the electron power partition increases with the toroidal mode number, while the hydrogen power partition decreases. This is due to stronger electron absorption with increasing toroidal mode number [9,10,45].
The ion absorption profiles for different toroidal mode numbers are shown in figure 9(a), where the results from FEMIC1D (solid lines) and TOMCAT (dashed lines) are compared. Good qualitative agreement was found with respect to the magnitude and absorption location, with some minor variations caused by different interference pattern. Figure 9(b) compares the electron absorption for n f =20. Two important observations can be made. First, the difference in the interference patterns is also observed in the electron absorption profiles. Second, TOMCAT predicts higher electron absorption in the interval   R 2.0 2.3 m. The difference is due to the IBW damping. In FEMIC1D, the IBW is completely absorbed at R=2.2 m. In contrast, the weaker damping in TOMCAT allows the IBW to propagate to the plasma edge at R=2.0 m, which increases the total electron damping in the interval 2.0<R<2.3 m. A similar behavior   was observed for the other toroidal mode numbers, which are not shown here for brevity. The electron absorption is described by TTMP and electron Landau damping. For the fast wave, these two absorption mechanisms are coherent and connected by the susceptibility elementsĉ yz e , andĉ zy e , . This is illustrated in figure 10 for the case with n f =20. Both the TTMP and electron Landau damping are positive definite, where the TTMP is a factor of two larger. The contribution of the cross terms is negative definite and roughly cancels the contribution of the TTMP, which leaves a net absorption similar to the electron Landau damping [1,45]. Figure 11 shows the absorption profiles of the IBW for different toroidal mode numbers, as predicted by FEMIC1D. The location of the maximum absorption is dependent on the    toroidal mode number. The lower the toroidal mode number, the further the IBW propagates before being damped. Similar conclusions can be drawn by studying the kinetic flux of the IBW in figure 12, where we also observe that the kinetic flux decreases as the wave propagates in the negative major radius direction. In these calculations, the total electron absorption of the IBW amounted to 0.63%, 0.44% and 0.21% of the total absorbed power for the toroidal mode numbers 15, 20 and 27, respectively. Thus, the fraction of mode conversion increases for lower toroidal mode numbers.
Since the IBW is an electrostatic wave, its Poynting flux is negligible. Therefore, the local absorption P abs is balanced by the divergence of the kinetic flux ·  T. This is illustrated in figure 13, which shows that these two quantities balance each other in the interval   R 2.2 2.5 m.

The correction current
The algebraic susceptibility c j in equation (7) can describe the physics of the fast wave quite well (including some FLR effects of the fast wave). However, this susceptibility model cannot describe the physics related to non-local effects, in particular, the mode conversion and the propagation and damping of the IBW. These effects, which are completely missing in the algebraic susceptibility model, are described by the correction current dJ ind . Figure 14 shows the components of the correction current. The oscillation and amplitude variations in δJ x and δJ y correlates with the oscillations of the IBW. The wavelet spectrogram of the δJ x component is shown in figure 15(a), which is useful in order to further understand the behavior of the correction current (the wavelet spectrogram for δJ y shows similar information and is not shown for brevity). In the spectrogram, we can see that the wavelet bases that contribute     to the induced current density coincide quite well with the wavelet bases of the IBW in figure 6. Another interesting point is that the wavelet bases related to the fast wave branch range from very small to negligible. This indicates that the algebraic susceptibility model is a good approximation to describe the propagation and damping of the fast wave, as the correction current does not add any significant corrections to the fast wave. The parallel current density component δJ z is shown in figure 14(c), and its wavelet spectrogram in figure 15(b). This component adds corrections to both the fast wave and the IBW. The corrections for the IBW are related to theĉ zz j , component and are rather small. This is expected, since the dependence of k R in theĉ zz j , is very weak. In contrast, large amplitude corrections for the fast wave are observed. These corrections are also caused by FLR effects, and affects the fast wave polarization of the wave in the yz plane. The effect originates from the susceptibility elementˆĉ c =zy j yz j , , , where the fast wave polarization is described bŷĉ c » -E E z y zy j zz j , , [45] (a similar effect is caused by thê c zx j , term, but is negligible here). Sinceĉ zy j , is proportional to the wave number k x , the response is dependent on the sign of the wave number, which causes the effect. However, this effect does not modify the power absorption, since the cross term in the power kernel,ĉ E E z zy j A y , , becomes an even function of k x (the superscript A denotes the anti-Hermitian part). Thus, calculations that do not account for this FLR effect may still obtain a reasonable power balance.

Convergence
The residual norm is based on the electric field components in the Anderson Acceleration routine, which is given by Convergence is achieved when the residual norm f is smaller than the absolute tolerance, which is set to 10 −3 . The evolution of the convergence is shown in figure 16. Overall, as long as the FE and wavelet grids have good resolutions, the Anderson Acceleration fixed point iteration scheme exhibits good convergence properties. However, the number of iterations to achieve convergence vary significantly. It was found that the number of iterations required for convergence are related to the number of oscillations in the correction current. Hence, the longer distance the IBW propagates, the more oscillations will be required in the correction current, which increases the number of iterations.

Discussion
In general, the propagation and damping of the fast wave as predicted by the FEMIC1D code is consistent with the predictions of the TOMCAT code. The observed difference is the damping of the IBW, where FEMIC1D predicts stronger electron Landau damping compared with TOMCAT. The reason behind this is that the FEMIC1D code has an all order FLR description of the dielectric tensor, while the TOMCAT code uses a truncated FLR expansion. A truncated FLR expansion works well when the Larmor radius of the hydrogen population is small, i.e.  l 1 H . This model works well for the fast wave and for modeling mode conversion. However, as the wave propagates away from the mode conversion  zone, the wavelength of the IBW becomes shorter and at some point it will become comparable to the ion Larmor radius, so that l~1 H . The consequence when the smallness assumption is broken is that the electron Landau damping becomes significantly weaker compared with an all-order FLR model, as was discussed in [41].
To evaluate equation (6) and account for the non-local response, a wavelet approach was used. The evaluation can also be done with Fourier bases, however, the Morlet wavelet approach has two advantages over the Fourier approach. The first is that Morlet wavelets are localized in space while having a narrow wave number content. Morlet wavelets are therefore suitable for describing waves in inhomogeneous media. The second advantage is that the wavelet transform does not require a periodic function on a finite grid. Fourier decomposed functions are required to be periodic in order to avoid the Gibbs phenomenon. Instead, wavelets have the problem of producing edge artifacts, which occur when Morlet wavelets are close to the boundaries. In this work, the issue is solved with the matched layer technique (see section 3.5), which is computationally inexpensive.
The main challenge of the iterative wavelet FE scheme is the number of iterations required to achieve convergence, which makes the numerical method slow. The number of iterations appears to correlate with the number of wavelengths in the correction current δJ ind . In [27], a wavelet spectral method was used to solve the wave equation and study fast wave heating in a toroidal plasma, which could account for non-local response in the parallel dispersion (i.e. the up-and downshift in the parallel wave number). The wave equation was solved by describing all differential operators in wavelet space, including the curl-curl operator, i.e. solving the implicit equation where c 0 is the cold susceptibility tensor. The required number of iterations to achieve convergence was ∼100. In [24], a similar wave equation was solved, but using the iterative wavelet FE scheme. For this case, the required number of iterations to achieve convergence was reduced to ∼10. Hence, the slow convergence in the present model is related to the fact that the wave equation describing the propagation of the IBW is treated using wavelets. Therefore, a possibility to reduce the number of iterations is to move the lowest order FLR effects from the correction current using operator splitting and include them in the FE method, and let the correction current contain higher order FLR effects. In FEMIC1D, the most time consuming part is the calculation of the inverse continuous wavelet transform of the   (44)) evolution with the iteration number for the three toroidal mode numbers. correction current (equation (40)). One approach to improve the performance of the inverse wavelet transform is to remove the redundancy in the wavelet representation of the functions, in particular, for large wavelets. For example, the discretization in the shift parameter b can be improved by using a coarser resolution for larger wavelets, while maintaining fine resolution for smaller wavelets only.

Conclusion
In summary, we have developed a new numerical method based on an iterative wavelet FE scheme, which has the capability to solve the wave equation and take into account FLR effects to all orders. To verify the method, we studied a case of onedimensional fast wave heating at the second harmonic ion cyclotron resonance of hydrogen and study mode conversion to IBWs. Comparison with simulations performed with the TOMCAT code shows good agreement in regions where the Landau damping of the IBW is not significant. The numerical method can also account for FLR effects related to the polarization of the fast wave in the yz plane, which are caused by theĉ zy j , term in the susceptibility tensor. The main advantage of the proposed numerical scheme is that it relies mainly on the FE method, which is suitable for solving the wave equation in complex geometries. In particular, the method can use the FE method in the low-density and antenna regions in a tokamak, where spatial dispersion effects are either weak or absent. Iterative addition of spatial dispersion effects are applied locally (i.e. the plasma domain), since the correction current δJ ind is evaluated on the wavelet grid, which is smaller than the FE grid. Another advantage of the method is that the problem solved with the FE method does not have to be a rigorous description of the plasma response, as iterations are used to approach the rigorous description.
The main challenge of the iterative wavelet FE scheme is to reduce the number of iterations required to achieve convergence (which is left for future work). One way to achieve this is to move the lowest order FLR terms from the correction current and treat them using the FE method, and only iterate on higher order FLR terms.