3D Full-Wave modelling and EC mode conversion in realistic plasmas

The wave physics of O-X conversion in overdense W7-X plasma is discussed. For this study, a new 3D, cold plasma full-wave code has been developed. The code takes advantage of massive parallel computations with Graphics Processing Units (GPU), which allows for up to 100 times faster calculations than on a single-CPU. A 3D calculation of the O-X conversion is demonstrated. We discuss limitations of the mode conversion scenario within the capabilities of the existing ECRH system in W7-X, and demonstrate an optimised conversion scenario in which the launching antenna location is altered. The conversion eﬃciency of the optimised scenario is predicted to be >85%.


Introduction
Electron Cyclotron Resonance Heating (ECRH) is the main plasma heating mechanism in Wendelstein 7-X (W7-X) Stellarator. It is provided by 10 gyrotrons at 140 GHz (corresponding to the second harmonic cyclotron resonance at 2.5 T) with the power of 1 MW each. X-and the O-modes were successfully used in a wide range of operation scenarios: X-mode for low and moderate densities (up to the cutoff at 1.2 · 10 20 m −3 ), and O2-mode for higher densities (up to 2 · 10 20 m −3 ). Possible operation at yet higher densities would involve double modeconversion from O-to slow X-and to Bernstein-mode, i.e. an OXB-scenario. The physics of O-X conversion is outside of applicability of the routinely used geometrical optics approximation (WKB-theory) and should be considered within a full-wave approach.
This work reports on the development of a new 3D cold plasma full-wave code, named CUWA. The code utilizes the Finite Difference Time Domain (FDTD) technique [1,2] and has an interface with the ray-tracing code TRAVIS [3]. The computation domain is "minimized" around the WKB-trajectory obtained from the ray-tracing code by means of the Convolutional Perfectly Matched Layer (CPML) technique [4]; the background magnetic field is recovered from pre-computed 3D equilibrium data. The code takes advantage of massive parallel computations with Graphics Processing Units (GPUs), which allows for up to ×100 acceleration over a single-CPU.

Numerical model
FDTD (also known as Yee's method) is a standard technique for full-wave simulations in various branches of physics [1,2]. In fusion plasmas, it has been successfully applied to study wave propagation, O-X conversion [5], reflectometry [6] and the scattering of EC waves on plasma * e-mail: pavel.aleynikov@ipp.mpg.de turbulence [7]. In this work we mostly follow the FDTD method as given in Chapter 11.3 of [2], with two important differences: first, we modify the current discretisation scheme and, second, we amend the differential operators of Maxwells equations in the boundary layers in accordance with the CPML technique [4] (see below).
The system of equations been solved is Maxwell's equations with addition of a "cold plasma" response current equation (J) given by the electron law of motion: where ω 2 p = n e e 2 /ε 0 m e is the plasma frequency, ν is electron collisional frequency and ω c = |e|B 0 /m e is the cyclotron frequency corresponding to the background magnetic field, B 0 .
In the FDTD method the field components are discretised on staggered grids in space and time. Note that is it common to discretize the plasma-response curren, J, in such a way that its J x , J y and J z components are co-located in space [2,5]. This facilitates the current update calculation. However, we find the stability of this scheme to be unsatisfactory, in particular when applied to the CMPL region. Therefore we have implemented a slightly more computationally demanding but more stable scheme where J is discretized in the same way as the wave electric field, E. This implies that the the current update equation will involve interpolation of its components.

Perfectly matched layer
Minimization of the computation domain plays a critical role in efficient 3D computations. A physical absorber layer with a non-zero collisionality (ν > 0) is often used in plasma physics computations to truncate the computation domain [5,7]. In such a case one would need to amend the computation domain with layers of "many-wavelengths" thickness at each boundary ensuring smallness of reflected signals. This may become problematic in 3D calculations when the domain of interest itself has the size of a few wavelengths. Amending such a domain with six boundary layers may significantly reduce the overall efficiency.
A so-called CPML boundary [4] is instead implemented in our full-wave code to truncate the computation domain. An efficient CPML requires only a few extra Yee cells in each direction -i.e. a fraction of a single "wavelength". CPML also requires two auxiliary (virtual) 3D fields (ψ) for each component of E and B, i.e. 12 auxiliary fields in total. However, these auxiliary fields are defined only within the narrow CPML layers, resulting in a negligible overhead in the practical setup.
The conceptual idea behind CPML [4] can be summarised in the following way (further details and a rather complete overview of various PMLs can be found in the book [8]). CPML is an efficient numerical implementation of Complex Frequency Shifted PML [9], in which the differential operators in Eqs. (1) are replaced with the convolution: where L −1 {F(ω)} denotes inverse Laplace transform and (κ u , α u , σ u ) are free parameters accounting for coordinate stretching in real and complex plains and an "effective conductivity". PML given by Eq.
(2) introduces no reflections because it can be viewed as a vacuum with complexstretched coordinates. In other words, with Eq. (2) the solution of Eq. (1) is analytically continued into the complex plain with a subsequent mapping of the complex coordinates back to real space but with a "complex material". In order to demonstrate CPML performance for the FDTD realisation described above we design a benchmark which represents a general anisotropic scenario. In this scenario, a divergent Gaussian beam propagating through the plasma escapes from the computation domain through a strongly nonuniform boundary. The entire computation domain is surrounded by the narrow CPML layer. This is shown in Fig. 1, the wave is launched vertically from the origin (0, 0). The plasma density ranges from 0 to approximately 80% of the O-mode cut-off with a gradient directed diagonally (bottom right to top left). It can be seen that at the top PML layer the wave-front incident angle varies from almost normal to quite oblique, ensuring extended benchmarking. The resulting reflection error is estimated as a maximum value of the electric field amplitude recorded in the bottom-right quarter of the domain. Note that the maximum relative error reported in the original work [4] is ≈ 5 · 10 −4 . Similar reflection is observed in our benchmark case (a wave with the cyan contours [10 −4 ; 10 −3 ]). Note the presence of another wave going out through the top right corner (contours [10 −5 ; 10 −4 ]). This spurious wave is associated with an inaccuracy of the Gaussian beam setup at the bottom boundary. Frequently, CPML outperforms the initial condition accuracy (for instance in a typical case where PML lies in vacuum) then the latter becomes the dominant noise source (i.e. introduces ∼ 10 −4 error in field amplitude). Note that the relative errors are squared when the intensity or mode power is of interest. The above example is calculated with the resolution of 12 Yee cells per vacuum wavelength and a CPML size of 12 cells.

Mode conversion examples
One of the applications of our code is the calculation of the O-to X-mode conversion efficiency in 3D plasmas. Figure 2 shows an example of a full-wave 2D calculation of the O-X conversion process in a realistic W7-X plasma. The O-mode, launched from vacuum (bottom right corner of the left plot), propagates into the plasma and is partly reflected off the corresponding cut-off layer, where a part of the O-mode is converted into a slow Xmode, which propagates toward the upper-hybrid resonance. Eventually, this X-mode will be converted to a Bernstein wave. The reflected part of the O-mode has a "hole" (due to the non-uniform O-X conversion). The perpendicular component of the refractive index is computed with the use of a windowed Fourier transform and shown in Fig. 2 (right). The incoming O-mode, the reflected O-mode, and the slow-X mode (which settles near the Upper Hybrid resonance layer) can be clearly distinguished in this diagram, which appears to be in a good agreement with expectations from WKB theory (i.e. dispersion curves), validating FDTD solution. Note that the slow-X mode is eventually dissipated numerically without leaving the computation domain. This dissipation is of purely numerical nature. It happens when the wavelength of the mode becomes shorter than 10 computational cells. Figure 3 shows an example of a 3D calculation similar to the above scenario (the wave is launched from the bottom of the domain in y direction, the slow X-mode is contained within the computation domain, the reflected Omode power distribution is shown on the plane x = const ). The decimal logarithm of the time-averaged E 2 is colourcoded. The sequence of dark-colored points represents the corresponding WKB trajectory for the reference ray. It is evident from the peculiar cross-section of the reflected wave that the 3D geometry needs to be appropriately accounted for.
To estimate the power conversion efficiency, we compute surface integrals of the time-averaged Poynting vector (i.e. the intensity) over every boundary of the computation domain. In the above case, the OX conversion efficiency can be calculated as the difference between the intensities of the input (S in ) and reflected (S ref ) waves. An estimation of the OX conversion efficiency is made as follows for W7-X-relevant ECRH beam parameters and plasma geometry. A Gaussian beam approximation is used as an initial condition. A beam of 2 cm radius and a phase front curvature of 1 m approximates the ECRH beam at the plasma entrance. The density profile is chosen arbitrarily (with a normalized gradient of ≈ 20m −1 near the cutoff). The antenna launcher location was fixed, corresponding to the upper port "A, C1" in W7-X nomenclature. Then, the launcher azimuth and altitude angles were varied to find direction for highest OX conversion. We find that for the given conditions (equilibrium, density profile, launcher position and beam focusing) the maximum conversion (S in − S ref )/S in reaches 48% (the middle tile in Fig. 3 (right)).

Optimised conversion scenario in W7-X
It appears that ∼50% is a typical maximimum conversion fraction for realistic W7-X beam parameters (focusing and width were optimized) and launcher positions. Refs. [10][11][12] suggest that conversion can be improved by a careful design of the beam wavefront. In this report we pursue an alternative optimisation strategy.
We take advantage of the fact that stellarator equilibrium "offers" a wide variety of magnetic field curvatures on a cutoff surface (flux surface), and we search for a location which would minimise the thickness of the evanescent layer, weighted over the beam radial profile. The thickness of the evanescent layer at the O-mode cutoff surface is estimated using the results of analytical theory [13], where the thickness is taken to be the distance (in wave-vector space) between the WKB solutions of the two modes. This estimate is known to agree very well with calculations in a flat 2D case [5].
The result of this calculation is demonstrated in Fig. 4  (left). The function representing the Gaussian-weighted gradient of the thickness of the evanescent layer is plotted on the cutoff surface, which is parametrized by the poloidal and the toroidal angles. A few minima of this function are marked with red dots. A different optimal value of parallel refractive index N opt || corresponds to each minimum due to the variation of the magnetic field strength on the surface.
It is easy to understand the "meaning" of these special points. We recall that the efficient conversion occurs when the "turning" points for the wave-packets of the two modes (O-and slow X-) are close to each other. Formally, the conversion condition is satisfied at the intersection of two surfaces: the O-mode cutoff surface and the X-mode turning surface parameterised by the plasma frequency (ω p ), the cyclotron frequency (ω c ), and the wave parallel refractive index N || . For a fixed equilibrium the in-  tersection curve (3D-shaped, in general) will vary with N || . Fig. 4 (middle) shows this intersection for N opt || ≈ 0.564, where the surfaces cross each other obliquely and thus the evanescent layer is narrow only close to the intersection (if the beam has a constant N || structure). When N || corresponds to one of the optimised locations, the intersection has a more complex structure (a saddle point in an example of Fig. 4 (right)), ensuring minimisation of the evanescent layer weighted over a finite Gaussian beam.
The location of the conversion region in plasmas (and thus N opt || ) uniquely defines the O-mode ray which arrives to this point within the WKB approximation. Therefore the beam launching location and angles can be computed with a backward ray-tracing calculation (staring from the conversion region). When the initial conditions for the beam outside of plasma are known, the full-wave calculations are done to predict the conversion efficiency. The maximised conversion efficiency predicted with our fullwave code for a case of Fig. 4 (right) reaches a remarkable 85%.
Note that such an algorithm for an identification of an optimised conversion scenario implies that the wave launcher location can be chosen freely. This is generally not the case in W7-X. Therefore our future plans include a study in which various W7-X equilibria are analysed in order to identify those where the optimal conversion regions are accessible within the capabilities of the existing ECRH system.