On all - FLR order numerical modelling of RF wave propagation and damping

. Determining the fate of waves in hot magnetized plasmas in magnetic confinement fusion machines close to arbitrary cyclotron harmonics retaining finite Larmor radius effects requires solving the relevant integrodifferential wave equation for waves of arbitrary wavelength. Anticipating exploitation in as realistic geometry as possible this requires massive computer resources. Work is ongoing to attempt reducing the required CPU memory and time. As Morlet wavelets are localised in x-space as well as k-space, they potentially offer a means to solve the relevant wave equation at reduced CPU requirement cost. Encouraging first results are obtained.


Introduction
Ion cyclotron resonance heating (ICRH) is a routinely used method to bring plasmas to fusion relevant temperatures in magnetic confinement fusion machines [1]. Properly modelling the dielectric response is a challenge both from the physics and the computational point of view. More often than not simplifications are made, first of all by truncating the dielectric tensor at the leading order terms in finite Larmor radius (FLR) corrections and secondly by simplifying the geometry. A new solver is under development to help understand the ICRH wave propagation and damping in multiple dimensions. It is inspired on the basic philosophy pioneered by Jaeger [2] for the AORSA code, and intends to provide a rigorous computation of the wave polarisation at arbitrary cyclotron harmonics. It will be useful both to study hot plasma wave propagation and damping effects in fusion relevant plasmas, and to add kinetic effects to the description contained in codes describing the electric field close to the antenna. The key idea at this initial exploratory stage is to assess the importance of higher order finite Larmor radius corrections while keeping the zero order distribution functions simple (Maxwellian) and the geometry prescribed.

Finite Larmor radius effects
A hot plasma admits an infinity of wave modes. Most of the key physics is captured by the "FLR2" models (retaining up to second order FLR kinetic terms i.e. corrections in k p ρ where k p is the perpendicular wavelength and ρ the Larmor radius) routinely used to describe the wave propagation and damping. Some physics aspects are outside this approach's scope, however. Examples are heating at higher cyclotron harmonics (N>2 where N is the cyclotron harmonic), wave absorption by very energetic populations and the fate of power carried by short wavelength waves (in particular the waves excited by the fast wave at the ionion hybrid layer).
The FLR2 model resting on the smallness of the k p ρ, upgrading the model to avoid errors arising from violations of this assumption necessitates solving the wave equation in its integro-differential form. Solving such an equation requires dedicated tools.
Jaeger solves the relevant equation directly in k-space, relying on the collocation method [2]. In his approach, the continuous integral in k-space allowing to account for the dielectric response of the hot plasma to the wave is replaced by a discrete set of Fourier integrals. In order to be rigorous, the considered k-spectrum should cover all physically important modes. Although proven to be extremely powerful, Fourier analysis has some drawbacks: since all Fourier modes are coupled, the system matrix of the corresponding linear problem is full, requiring vast computer resources. Another drawback is that Fourier analysis is ill-suited to describe evanescent modes since its building blocks (although forming a complete base and thus mathematically perfectly suited) represent propagative i.e. oscillating waves. It takes a lot of oscillatory base functions to approximate exponential growth or decay sufficiently accurately.
The ultimate goal of this work is to solve the hot plasma wave equation in 2 or 3 dimensions accounting for the realistic geometries of tokamaks and stellarators. Two roads are being pursued in parallel: (i) a "brute force" Fourier technique that puts the computational effort on supercomputers as it requires inversion of massive, full matrices, and (ii) a technique that explores the possibility of reducing the computational effort by semi-analytically solving some of the required integrals by hand and exploiting the fact that the net effect of very fast oscillations potentially offers a justification for reducing the full matrix to a sparse(r) one. Whereas the first is directed at exploring available numerical techniques to solve large linear systems swiftly and accurately, the present paper focusses on the reformulation of the equation in a suitable way attempting to avoid the need to invert a large, full matrix by transforming it in a sparse(r) one. Prior to exploring multi-dimensional application, tackling the integrodifferential equation in 1D is thought to be a first necessary step. low order polynomial. The equation at hand is solved in each subdomain for a complete set of "base solutions" consistent with the equation. Hence it represents all possible linearly independent solutions for the degree of accuracy imposed by the chosen base polynomials. The remaining degrees of freedom are exploited to impose the connection of the physical solution across interfaces. Although this technique was demonstrated to be very versatile and fairly easy to implement, it proved unable to provide converged solutions for the wave equation at hand.
A wavelet approach has now been attempted. Wavelets have the advantage that they are local, a characteristic they share with finite element base functions. In an elegant way, this potentially yields sufficiently sparse rather than full matrices for the non-local wave equation at hand. Moreover, certain classes of wavelets resemble wave packets and are therefore particularly suited to describe waves. In the magnetic confinement fusion domain, this technique was successfully pioneered by Vallejos et al. [5].

Wavelet approach
Morlet wavelets are a much adopted kind of wavelets. They are a family of functions and a suitable constant. For wave applications, it is more logical to use = 1/ rather than as a variable. This has the supplementary advantage that it removes a singular factor 1/a 2 from the integrand of the inverse wavelet transform, which is defined as for a function . The Morlet wavelet is localised in space around = and in k-space around / . An example is shown in Fig. 1. The top graph depicts the wavelet and the bottom graph its Fourier transform.
The method can be illustrated via a simple example: the solution of the wave equation with 2 boundary conditions of the type at the extremities of the interval of interest can be found by introducing the definition of the wave transform in the equation and the boundary conditions. The coefficients Ε( , ) are unknown while the base functions and their derivatives are. Similar to the technique adopted by Jaeger, the continuous Fourier integrals are approximated by discrete Fourier sums on a grid ! , ! and the equation is imposed at as many grid points as there are coefficients Ε ! , ! . At the extremities of the interval of interest, the boundary conditions are either imposed via 2 Lagrange multipliers or -simpler -substituted for the differential equation. The resulting linear system can be resolved by standard matrix inversion techniques. Since the wavelet base is localised, the system matrix is sparse: Introducing a suitable clipping criterium, contributions from 's too far from the reference position can safely be neglected. As an illustration, the solution of the Airy equation Ε′′ + Ε = 0 for wave incidence from the right is depicted (' = / ). Both regions of propagation and evanescence are treated to satisfaction.

Hot plasma wave equation
The homogenous plasma dielectric response of a hot plasma to an electromagnetic wave is well known [6,7].
In inhomogenous plasma the wave spectrum contains many coupled modes, each of which have a different response. It is customary to account for the inhomogeneity in an approximate way locally using a "quasi-homogeneous" assumption and writing the total RF current density as an integral composed out of contributions belonging to the various modes of the spectrum but assuming the dielectric response for a given k can simply be written as the response in a uniform plasma. The relevant equations are then given in the above given references. The hidden subtlety is the fact that the equation of motion is solved locally pretending the confining magnetic field is constant. Inside the plasma, the wave equation can then be written as in which ( ) is the electric field vector, ! its Fourier transform and , the homogeneous plasma dielectric tensor; k o =ω/c, in which ω is the angular frequency of the driver and c the speed of light. Thanks to the available expressions for this tensor, the formulation allows to study wave propagation and damping at arbitrary harmonics and for any wave type. The expressions for the dielectric tensor are known explicitly for bi-Maxwellian distributions with a parallel drift. Relying on numerically evaluated derivatives of the derivatives the corresponding expressions for the dielectric response are equally available for exploitation for non-Maxwellian distributions, as was e.g. demonstrated by recent work involving AORSA [8].
Substituting the Fourier transform of the electric field by its definition one obtains the non-local wave equation which can -be it at the price of more bookkeeping -be solved in the same way as described in the simple example. For 1D application in the -direction, the equation becomes in which = ! . A truncated Taylor series expansion of the dielectric tensor has been made in k= − and while x = ! − and ξ = − . Although the k and x integrals run from −∞ to +∞, the integration interval can be truncated in practice since the wavelets are localised. This is necessary to ensure the adopted Taylor series expansion is justified. The ! integrals can be integrated analytically. An example is shown in Fig. 3. The !!! coefficients are the coefficients obtained when rewriting ∇×∇× ( ) in terms of the derivatives of the field components.
To ensure the system matrix becomes sparse, negligible contributions to it should be omitted. It can easily be seen that if | | is small the contribution should be retained. But when it is large, looking at the shape of the base function and only retaining contributions larger than the inequality (−2 ) !/! < | | can be adopted as clipping requirement.

Numerical example
Although the testing of the method is still fully ongoing, a preliminary example is provided. The case examined is one of the key RF heating schemes in JET: fundamental cyclotron hydrogen minority heating in a deuterium plasma. The following parameters were taken: The magnetic field strength is B o =3.45T, which allows for core cyclotron heating when adopting f=51MHz as the driver frequency. Dipole phasing was considered but only the dominant component in the antenna spectrum (n tor =27) was retained. The central density is 3x10 19 /m 3 and the core temperature 3keV. The minority concentration is assumed to be 5%. The JET major radius is R o =2.97m and a minor radius of a p =0.9m was considered. The equation is solved in slab geometry but retaining the toroidal curvature; poloidal field effects are ignored. The poloidal (y) and toroidal (z) directions are assumed to be ignorable; k y =0 was chosen. The relevant fast wave dispersion roots is depicted in Fig. 4; x=R-R o . The fast wave becomes evanescent at the far high field side (left in the picture) and has a confluence with the Bernstein wave in the plasma centre. The dominant damping is expected close to the minority cyclotron layer (ω=Ω H at x=0.14m). The electric field components are depicted in Fig. 5. A fast wave is incident from the low field side. Since the Bernstein wave (excited near the ionion hybrid layer) is electrostatic, no trace of short wavelength behaviour is seen on the E y field component but it is visible in the E x component. The component of the electric field parallel to the static magnetic field is small in the whole interval. At the far high field side the wave is reflected from the cutoff across which its amplitude decreases exponentially. The presence of the Berstein wave can also be seen in the Fourier spectrum of the field; see Fig. 6. The biggest peak lies at k x =-25/m and corresponds to the fast wave incoming from the low field (propagating towards negative x); the slightly smaller peak with the opposite sign is the reflected fast wave. After being born at the mode conversion layer the Bernstein wave propagates leftward. This wave being a backward wave with a wavelength shorther than that of the fast wave, its signature is seen to the right of the fast wave root at k x >25/m.   Although the obtained solution contains many features that are in line with expectations, there are various clear hints that the solution is not yet converged. A closer look at the power balance reveals that the majority absorption dips slightly below zero, which is unphysical for a plasma in which all distributions are Maxwellian. Also the fact that the electrons do not take part at all in the absorption is highly unusual. Electron absorption scaling with k⊥ 2 when adding the electron Landau and transit time magnetic pumping contributions [9,10], the absence of electron absorption suggests that the considered spectrum not containing high enough k x -modes. More in depth analysis confirms this is indeed the case: see the inset for the corresponding evaluation using the TOMCAT differential equation solver [11] which shows a small but finite electron absorption exhibiting both long (fast wave) and short (Bernstein wave) wavelength structure.

Final remarks
In this paper the wavelet technique was adopted to solve the hot magnetized plasma integro-differential wave equation. As is done by Jaeger for AORSA, it maximally exploits the dielectric tensor expressions available in the literature. The approach proves promising but though more localised and thus yielding sparse matrices, the double integral on and is somewhat more CPU demanding than was anticipated and necessitates further optimization. Among the techniques under examination are the exploitation of pushing the analytical evaluation of the coefficients needed in the wave equation further and attempts to push down the CPU requirements by avoiding to invert the system matrix adopting e.g. iterative solvers for the linear system.