Benchmarking full-wave codes for studying the O-SX mode conversion in MAST Upgrade

. Three full-wave codes for simulating microwave propagation and O-SX mode conversion in magnetized plasma are described and compared. Their feasibility to investigate mode conversion processes and obtain conversion e ﬃ ciencies for parameters relevant for a potential MAST Upgrade 28GHz electron Bernstein wave heating scenarios is explored.


Introduction
Electromagnetic waves in the microwave range of frequencies play an indispensable role in fusion plasmas for heating [1] and diagnostic purposes [2]. One of their advantages are the small space requirements for in-vessel components compared to other heating sources which is of special importance in a potential reactor where the in-vessel wall will be covered with Tritium-breeding blankets [3].
In conventional tokamaks & stellarators the microwave energy is usually deposited where the angular frequency of the injected microwave ω 0 equals the electron cyclotron frequency ω ce = eB 0 /m e , with e the elementary charge, B 0 the background magnetic field, and m e the electron mass. If the electron plasma frequency ω pe = n e e 2 /( 0 m e ) (with n e the electron plasma density and 0 the vacuum permittivity) exceeds ω ce , the wave is in cut-off and an injected O-mode can no longer reach the electron cyclotron resonance (ECR) layer. Such a plasma is often referred to as overdense. One solution is to double the frequency of the injected microwave and heat at harmonics of the ECR. Second harmonic heating is a standard scheme in fusion plasmas but higher harmonics are rarely used as the absorption decreases with increasing harmonics number.
Compared to conventional tokamaks, spherical tokamaks exhibit a higher plasma-β [4], resulting in highly overdense plasmas. Heating at harmonics becomes ineffective due to the high harmonics number. A solution is provided by the electrostatic electron Bernstein wave (EBW) which has no high-density cut-off, is very well absorbed at the ECR even at high harmonics and can very efficiently drive toroidal net currents [5]. This makes them an attractive candidate for spherical tokamaks, where noninductive current drive techniques are mandatory due to a small (or completely absent) central solenoid. * e-mail: koehn@igvp.uni-stuttgart.de The paper gives an overview of a set of numerical tools to be used to explore the feasibility of a MAST Upgrade 28 GHz EBW heating scenario. The numerical set-up used in the simulations is described in Section 2, including a brief description of the different numerical codes. Section 3 contains the benchmark studies, Section 4 investigates the cases of very steep plasma density profiles. A brief summary concludes the paper.

The simulation domain
A Gaussian beam with a frequency of f 0 = 28 GHz in O-mode polarization is considered in the simulations. Following e.g. Ref. [2], its electric field distribution reads with x the radial distance to the beam axis, z the axial distance to the beam waist, w 0 the size of the beam waist (here w 0 = 4 λ 0 ), R the radius of curvature of the wavefront, and φ 0 the Gouy phase shift. Note that w, R and φ 0 are all functions of z. A linear plasma density profile is used with a normalized density scale length of k 0 L n = 25, corresponding to with z n,start = 0.15 m. The background magnetic field is set to a constant magnitude of |B 0 | = 0.85 T, pointing into the y direction: B 0 = (0, B 0 , 0). Note that z is understood here as the direction into which the microwave beam is primarily injected (which looses of course its meaning when considering O-SX mode conversion which requires an oblique injection with respect to the background magnetic field).
magnetized plasma 0 θ opt Table 1. Scenario configurations in the toroidal plane, z 0 corresponds to the position of the beam waist, θ refers to the injection angle with respect to the y-axis, where θ opt is the angle resulting in maximum mode conversion efficiency, Eq. (10).
The simulation box is set to a size of N y = 80 λ 0 and N z = 30 λ 0 . The different scenario configurations are summarized in Table 1: they increase in complexity with scenario #6 finally corresponding to a case of optimized O-SX mode conversion. Note that those scenarios all treat propagation in a "toroidal plane", i.e. y corresponds to the toroidal direction assuming a conventional tokamak whose magnetic field points primarily into the toroidal direction.
To compare results between the different numerical codes, we will use cross sections of the squared wave electric field summed up over a few wave periods (i.e. time-averaged) along the y direction. Cross sections at z = {0, 5, 10, 15, 20, 25, 30} cm will be considered. For the cases where O-SX mode conversion takes place, we will also compare the conversion efficiencies.

IPF-FDMC
IPF-FDMC is a 2.5D finite-difference time-domain (FDTD) code [6] solving Maxwell's equations and the fluid equation of motion of the electrons on a 2D Cartesian grid. The code is written in C and uses OpenMP for parallelization on shared memory systems. It allows to simulate propagation of electromagnetic waves in a cold magnetized plasma. The evolution equations considered for the magnetic field B, the electric field E, and the current density J of the wave in a plasma equilibrium are: with c the speed of light andB 0 the unit vector into the direction of the background magnetic field. An electron collision frequency ν e is included in Eq. (5) as a dissipation mechanism. The code has been successfully benchmarked against cold plasma theory [7] and against the full-wave code EMIT-3D [8]. For the implementation of Eqs. (3)-(5), the standard FDTD scheme has to be complemented by a discretization scheme for the current equation (5). IPF-FDMC uses a "straightforward" way: first advance J x using the old values of J y and J z (in the cross product with B 0 ), then advance J y using the updated value of J x and the old value of J z , and finally advance J z using the updated values of J x and J y . As has been shown in previous works, this method is completely sufficient for a rather large set of problems [9,10]. More advanced methods exist for situations with extreme density gradients [11,12]. A Gaussian beam is injected into the grid as a soft source [6] with the electric field amplitude, as defined in Eq. (1), specified in an antenna plane. The beam waist w 0 and the axial distance z to the beam waist are input parameters. According to e.g. Ref. [2], w and R in the antenna plane are given in vacuum by the following equations: with z R = πw 2 0 /λ 0 the Rayleigh range. Some algebra yields the expressions for the antenna input used in IPF-FDMC: The simulation domain considered here is surrounded by non-radiating boundaries. IPF-FDMC uses a physical absorbing boundary layer which damps the wave electric field by multiplying it with a damping function d(x) whose strength depends on the position x in the absorber: with λ grid the wavelength in grid points, set to 64 grid points per vacuum wavelength, d absorb = 3 λ 0 the thickness of the absorbing layer. This method also yields the amount of wave energy damped in the absorber and thus indirectly the O-SX conversion efficiency as the converted SX-wave is dissipated in the vicinity of the upper-hybrid resonance (UHR) and not in the absorber.

EMIT-2D
EMIT-2D is a 2D full-wave, cold plasma code parallelised using OpenMP. It is a novel code, developed from EMIT-3D [8] using the FDTD method to solve Maxwell's Eqs. (3) and (4). The plasma response is accounted for via Eq. (5) (with ν = 0). Equations (3) -(5) are solved on a Yee grid with J being defined on the same spatial points as E. A leap frog scheme is used, i.e. E and B are calculated on alternate time-steps (with J being calculated on the same time-step as B) using the previous value of the other. As EMIT-2D is a 2D code, derivatives in the third dimension (chosen as the y dimension arbitrarily) vanish. Boundary conditions are required to solve the update equations on the boundary grid points. In EMIT-2D, the electric field at the boundary is set at zero. In order to avoid reflections from the boundary grid points, an absorbing boundary layer is implemented. The boundary layer is three vacuum wavelengths thick, and in it all components of the electric field are multiplied by a cubic function: where τ is the wave period divided by the timestep ∆t, d bound is the boundary thickness and d < d bound is the distance from the edge of the simulation domain. This ensures that power returned into the simulation domain is diminished by a factor of less than 10 −8 [13]. The spatial resolution is set to 50 Yee cells per vacuum wavelength. A Gaussian beam is launched to propagate in the zdirection using a soft antenna, meaning the electric field is added to whatever is already present. The Gaussian beam launched by the antenna is described by Eq. (1) multiplied with a time-harmonic variation, sin(ωt).

Fourier Full-wave simulation code
The Fourier Full-Wave (FFW) simulation code solves Maxwell's equations coupled with the electron fluid equations, using a Fourier method in the plane normal to the density gradient. Separation of scales in time of the unknown variables is used leading to a system governing the slowly varying in time complex envelope of the wave electric field. Such separation of time scales are frequently used when studying quasi-monochromatic waves with slowly varying amplitudes and phases, leading to, e.g., the non-linear Schrödinger equation in non-linear optics and plasmas, and to the Zakharov equation when studying the non-linear coupling to ion waves [14][15][16][17]. Since we are studying the linear behaviour of microwaves, the non-linearity has been "switched off" keeping the ions immobile. Exploiting that the amplitude and phase are slowly varying, it is possible to take much longer timesteps in time than the oscillation period of the microwave.
An implicit trapezoidal (Crank-Nicholson) scheme advances the solution in time, with the advantage of being unconditionally stable with respect to the time-step and second-order accurate. In the simulations performed, one time-step corresponds to approximately 60 oscillation periods of the microwave, and for our purposes it suffices to take 50 -100 time-steps to obtain the desired solution.
By Fourier transforming the system of equations in the plane perpendicular to the density gradient, the system separates to one 1D system per Fourier component. Simple out-flow boundary conditions are used at the antenna plane z = 0, where waves with the desired polarization are injected [16]. The remaining spatial variable parallel to the density gradient is resolved by 40 grid points per vacuum wavelength. To resolve the small-scale electrostatic waves resulting from mode conversion and resonant absorption of the EM wave a region around the resonant layer is resolved by 2000 grid points per wavelength. After solving each of the 1D systems in space and time, the resulting solutions are superimposed with appropriate weights and phase factors to obtain the full multi-dimensional solution.

Scenario 1
A vacuum, no plasma and B 0 = 0, is considered in this scenario. A Gaussian beam, see Section 2.1, is injected from the bottom, with the waist at the antenna plane, propagating parallel to the z-axis. To compare the results of the full-wave codes, cross sections along various z-positions are taken, see Sec. 2. These signals will be referred to as detector antenna signals in the following. As shown in Fig. 1, one can see that with increasing distance from the antenna plane at z = 0, the amplitude of the beam cross section decreases and the width increases (as expected).
For a more quantitative comparison, the detector antenna signals of the different codes are subtracted from each other. This requires to map the results of one code to the spatial grid of another code. As an interpolation technique, Piecewise Cubic Hermite Interpolating Poynomial (a 1D monotonic cubic interpolation) from Python's SciPy module is used [18]. The resulting differences are shown in Fig. 2. The differences between IPF-FDMC and EMIT-2D are very small which is not surprising as both codes are based on the same numerical method. A further reduction can be expected if the spatial resolution of one of the codes were to be increased. The difference between FFW and the other codes is larger by an order of magnitude which is, however, still small (and thus acceptable).

Scenario 2
In scenario #2, the position of the beam waist is shifted upwards to z = 0.2 m, resulting in a converging beam. The quantitative comparison is plotted in Fig. 3: the differences between IPF-FDMC and EMIT-2D are again very small. The differences between FFW and the two FDTD codes are again somewhat larger but still small.

Scenario 3
In scenario #3, the waist is located at z = 0 and the beam is injected at an angle of θ = 30 • with respect to the z-axis. The quantitative comparison of the detector antenna signals gives a similar picture as for scenario #1: as shown in Fig. 4, the differences between IPF-FDMC and EMIT-2D are very small and can be expected to be further reduced by increasing the spatial resolution. Although the difference between FFW and the two FDTD-codes is larger by an order of magnitude, it is still small and, looking at the spatial distribution of the differences, seems to be mostly due a small shift (or rather an offset) of the spatial grid of FFW with respect to EMIT-2D and IPF-FDMC.

Scenario 4
In this scenario we have introduced the electron plasma density profile as described by Eq. (2), but still have B 0 = 0. A microwave beam with the waist located at the antenna plane is emitted at an angle of θ = 30 • . Since there is no magnetic field here, no O-SX conversion can happen (for B 0 = 0, the two solutions of the dispersion relation degenerate and O-and X-mode are indistinguishable). There is, however, coupling to plasma oscillations, also known as Langmuir waves, possible if the wave electric field has a component along the density gradient [19]. Figure 5 shows a snapshot of the wave electric field obtained from the full-wave simulations. Small scale oscillations can be seen at the plasma cut-off layer: these are the aforementioned plasma oscillations. In the codes, such small-scale oscillations are eventually dissipated via a collisional damping term. Note that for the case of FFW there seems to be a second wave injected at a position of roughly y = 0.5 m. This corresponds to reflections at the bottom absorber from the "main" wave leaving the grid at that position. While this might look like an issue, it is only visible due to the logarithmic color scale and is negligible when calculating power absorption in the grid. It is, nevertheless, planned to further improve the bottom absorber to suppress these spurious reflections.
FFW yields an absorption at the plasma frequency layer of approximately 4 %, IPF-FDMC of approximately 2.1 %. A formula by Mjølhus [20] allows to estimate this value to 5 % considering a plane wave description and thus serving as an upper boundary which would be approached if beams with increasing waist would be used in our simulations. The difference between FFW and IPF-FDMC is due to the different implementation of the beam, see the corresponding discussion in Sec. 3.6.
The differences between the detector antenna signals of the codes are shown in Fig. 6 (note that there is no signal for z ≥ 0.2 m as the wave is in cut-off at the corresponding densities). While for most detector antenna positions the difference is small, it is more pronounced when being close to the cut-off. The reason is probably a small mis- alignment of the numerical grids along the z-axis. When being very close to the cut-off, the shape of the horizontal cross sections strongly depends on that position.

Scenario 5
In this scenario, a homogeneous background magnetic field B 0 = 0.85 T is introduced pointing into the ydirection. The microwave beam is injected at an angle of 30 • with respect to the z-axis. A pure O-mode is injected requiring an elliptically polarized beam [21]. Despite being far off the optimum injection angle for O-SX mode conversion, conversion is actually found in the simulations: a value of 1.1 % is obtained from IPF-FDMC, 1.9 % from EMIT-2D, and FFW yields 0.2 %, which gives an idea about the error-bars when comparing the codes.
Analyzing the differences between the detector antenna signals, see Fig. 7, one sees small differences between IPF-FDMC and EMIT-2D. The differences between FFW and the two FDTD codes are larger by a factor 2 but looking at their spatial distribution it seems to be due to a small misalignment of their spatial grids.

Scenario 6
The microwave beam is now injected at the optimum angle θ opt for achieving maximum mode conversion efficiency, corresponding to θ opt ≈ 43 • . A snapshot of the resulting wave electric fields is shown in Fig. 8. An enhancement in the vicinity of the UHR can be clearly seen. Another important feature is the "hole" in the reflected, non-converted part of the beam: the SX-mode content is "filtered" by the UHR (eventually dissipated via collisional damping in the codes), it cannot pass it. Only the tails of the Gaussian beams, which have non-optimum k-components, can pass the UHR as they cannot couple to the SX-mode. An O-SX conversion efficiency of 86.4 % is obtained from IPF-FDMC, 86.7 % from EMIT-2D, and 88 % from FFW. The higher value obtained from FFW is due to its implementation of a Gaussian beam: it is represented by a spectrum of k components where each component is in pure O-mode polarization. In the two other codes the Gaussian beam is defined by appropriate spatial wave electric field distributions with the tails having non-optimum polarization. Conversion efficiencies from FFW are therefore to be understood as an upper bound. The difference between the detector antenna signals, shown in Fig. 9, reveal a noticeable mismatch only in the vicinity of the O-SX conversion layer. This is likely due to the different dissipation mechanisms in the codes. Note that the conversion efficiencies agree, however, very well.

Steep density profiles
Plasma scenarios in MAST Upgrade relevant for a potential 28 GHz EBW heating system include L-and Hmode scenarios and different antenna positions. The resulting value of k 0 L n defines the sensitivity of the O-SX conversion on angular mismatch: higher values of k 0 L n (corresponding to more shallow profiles) are more sensitive to angular mismatch. An estimation based on a high elongation, low plasma-β scenario yields a range of k 0 L n = 2 . . . 25. Figure 10 shows the O-SX conversion efficiencies obtained from full-wave simulations for this range, showing good agreement between the different codes. Furthermore, one can clearly see a peak in the conversion efficiency at around k 0 L n ≈ 7. For smaller values (thus steeper profiles), the efficiency starts to decrease as now the SX-mode can tunnel through the evanescent layer defined by the UHR and the right-hand cut-off and couple to the FX-mode which simply leaves the plasma.

Summary
Three full-wave codes have been benchmarked against each other, where good agreement was found. All were shown to be able to simulate the O-SX mode conversion yielding very similar results. This makes them suitable tools for a feasibility study of an EBW heating scenario in MAST Upgrade, to be published in a subsequent paper. It was also shown that for too steep plasma density profiles, the O-SX conversion efficiency starts to decrease which might pose a threat to in-vessel components.

Acknowledgments
The authors are indebted to the efforts of the open-source software community.