Development and application of a ray-tracing code integrating with 3D equilibrium mapping in LHD ECH experiments

The central electron temperature has successfully reached up to 7.5 keV in large helical device (LHD) plasmas with a central high-ion temperature of 5 keV and a central electron density of 1.3×1019 ?> m−3. This result was obtained by heating with a newly-installed 154 GHz gyrotron and also the optimisation of injection geometry in electron cyclotron heating (ECH). The optimisation was carried out by using the ray-tracing code ‘LHDGauss’, which was upgraded to include the rapid post-processing three-dimensional (3D) equilibrium mapping obtained from experiments. For ray-tracing calculations, LHDGauss can automatically read the relevant data registered in the LHD database after a discharge, such as ECH injection settings (e.g. Gaussian beam parameters, target positions, polarisation and ECH power) and Thomson scattering diagnostic data along with the 3D equilibrium mapping data. The equilibrium map of the electron density and temperature profiles are then extrapolated into the region outside the last closed flux surface. Mode purity, or the ratio between the ordinary mode and the extraordinary mode, is obtained by calculating the 1D full-wave equation along the direction of the rays from the antenna to the absorption target point. Using the virtual magnetic flux surfaces, the effects of the modelled density profiles and the magnetic shear at the peripheral region with a given polarisation are taken into account. Power deposition profiles calculated for each Thomson scattering measurement timing are registered in the LHD database. The adjustment of the injection settings for the desired deposition profile from the feedback provided on a shot-by-shot basis resulted in an effective experimental procedure.


Introduction
Optimal injection settings of electron cyclotron heating (ECH) are essential for achieving the desired power deposition and reducing the stray radiation level in the vessel. In particular, realising on-axis heating can result in a high central electron temperature. Since plasma parameters, such as the profiles of electron temperature T e and electron density n e , vary during a plasma discharge [1,2], optimal injection settings have to be determined and adjusted according to those plasma parameters.
A ray-tracing calculation for the EC wave can predict the power deposition profile from the combination of the magnetic field configuration, the ECH injection settings, and the T e and n e profiles in a plasma equilibrium. A realistic calculation can be performed by taking detailed information into account. In particular, the refraction of rays in the peripheral region of helical plasmas such as in the large helical device (LHD) [3], which is the largest heliotron-type superconducting device in the world with the capability of generating numerous 3D equilibria in a variety of heliotron configurations, cannot be neglected due to a finite density gradient. In addition, since strong magnetic shear exists in the region, the polarisation of the injected millimetre waves has to be carefully adjusted in order to excite pure ordinary (O) mode or extraordinary (X) mode at the resonance layer. For an effective absorption of high power EC beams, the effects of refraction and pure mode excitation are also important in tokamak plasmas such as ITER or future fusion plasmas.
The precise evaluation of the deposition profile is also essential for the study of transport. For instance, modulation ECH (MECH) has been used to generate heat pulses, which are used to evaluate cross-field electron thermal transport [4][5][6]. In such studies, the precise identification of the deposition profile is quite important.
Recently, a rapid post-processing magnetic coordinate mapping system, or a three-dimensional (3D) equilibrium mapping system, has been developed and applied to LHD experiments [7]. The rapid post-processing equilibrium mapping is based on the selection of a best fit profile among thousands of equilibria pre-calculated with the VMEC code [8]. The Thomson scattering diagnostic system provides the input data for the mapping [9]. As a result, non-axisymmetric 3D equilibrium for every Thomson scattering time slice of every shot is offered rapidly after the shot during LHD experiments. The mapping data can be available after the discharge, i.e. the post-shot analyses to catch up with the experimental sequence of shots. This mapping system is now incorporated into the 'AutoAna' (auto-analysed data registration to the LHD analysed data server) system [10], which automatically processes physical raw data obtained from LHD experiments and registers analysed data such as rapid post-processing 3D equilibrium mapping of all the time slices of the Thomson scattering data during experiments.
Until now, conventional ray-tracing calculations using the ray-tracing code 'LHDGauss' without the rapid post-processing 3D equilibrium mapping system have been conducted for certain shots to study ECH power deposition [11][12][13].
Therefore, feedback of the injection settings to meet the required deposition profile was not available on a shot-by-shot basis during experiments, which often resulted in performing the experiments under an unintended heating condition (e.g. off-axis heating).
This paper reports the upgrade and the outcome of the raytracing code LHDGauss, which now includes the rapid postprocessing 3D equilibrium mapping system as well as the mode purity analyses for calculating the absolute values of ECH power deposition. Using this upgraded LHDGauss, we succeeded in increasing the central electron temperature T e0 of a high-ion-temperature T i plasma, and the high T e -T i parameter regime of the LHD plasmas was expanded.
This paper is organised as follows. Section 2 together with appendices A and B describe an overview of LHDGauss. In section 3, the calculation results and some experimental results are compared. In section 4, the successful results of increased T e with feedback-tuned ECH injection settings are presented. Section 5 summarises this paper.

Ray-tracing code 'LHDGauss'
LHDGauss is one of the multi-ray-tracing codes [14][15][16] that calculate ray trajectories by solving the eikonal equation under the WKB approximation and calculating the power absorption along the ray propagation. This code is used to calculate the propagation of EC waves from the injection antenna as a Gaussian beam in the LHD vacuum region into the LHD plasma and also to calculate the absorption of the waves in the plasma to obtain the power deposition profile. This code is currently dedicated to ECH experiments in the LHD. The ECH system in the LHD is summarised, for example, in references [17][18][19]. LHDGauss is unique in that the multi-rays are distributed so as to simulate the Gaussian beam adopted in the LHD ECH system, in that the rapid post-processing magnetic coordinate mapping system is adopted, and in that the mode purity, which is the ratio between the ordinary (O) mode and the extraordinary (X) mode, is determined by solving the 1D full-wave equation from the injection antenna to the absorption target point, where the effects of the n e profile and the magnetic shear at the peripheral region of an LHD plasma under a given polarisation are taken into account by using the extrapolated virtual magnetic flux surfaces. The effect of the refraction of rays at the peripheral region outside the last closed flux surface (LCFS) with a finite density gradient is also included in this code by use of the modelled density profile with the extrapolated 3D equilibrium mapping.
The outline of the calculation procedure of LHDGauss is described in figure 1.
First, the input files for the ray-tracing calculations are prepared in the calculation server for LHDGauss. The processed Thomson scattering diagnostic data, i.e. 'TSMAP', 'TSMESH' and 'TSWPE' are read from the LHD database server [7]. The rapid post-processing equilibrium mapping, based on VMEC equilibria, is applied for the equilibrium within the LCFS in the LHDGauss code. TSMAP is equilibrium mapping that relates the real coordinate R to the effective minor radius r eff using the T e and n e profiles. TSMESH also relates the real cylindrical coordinates ( ) Φ R Z , , to r eff for the LHD. The cylindrical coordinates are then transferred to the Cartesian coordinates (X, Y, Z ) for ray-tracing calculations. TSWPE fits the T e and n e profiles inside the LCFS with polynomial functions in r eff . A Gaussian function is added to a polynomial function to express an electron internal transport barrier with a peaked profile in T e [21,22]. Since the ( ) n r e eff profile is given not only inside the LCFS but also in the peripheral region (outside of the LCFS) to some extent where the Thomson scattering diagnostics are available, the ( ) n r e eff data outside the LCFS are fitted by exponentially-decaying functions connected continuously with the fitted functions inside the LCFS. Then, the 3D equilibrium map is virtually extrapolated far outside of the peripheral region to the first wall of the LHD by using an extrapolation method, which is summarised in appendix A. Although the n e profile in this region may not be a flux function, the extrapolated profiles are chosen to be consistent with the Thomson scattering diagnostics available outside the LCFS. It is recognised that the pressure on the magnetic flux surface can be expressed by the flux function in the LHD irrespective of the neighbourhood of the diagnostics position or EC-beam launch ports. A possible error regarding the extrapolation process arises from the reduced accuracy of the extrapolated points far from the LCFS due to lack of the points inside the LCFS. However, the error does not affect either the polarisation or the refraction of rays since the electron density in these regions are too low. Thus, the extrapolated 3D equilibrium mapping with the interface region between the plasma region and the vacuum region enables the modelling of the propagation of the EC waves in the entire region from the injection antenna to the plasma core.
All vacuum magnetic field components on given mesh points for most of the combinations of each magnetic field coil are pre-registered in the LHDGauss calculation server, and one of them is referred to before the ray-tracing calculation. Since the 3D equilibrium mapping does not provide the magnetic field information outside the LCFS, the vacuum magnetic field components are used for the ray-tracing calculation at present. The volume-averaged beta for the shots Figure 1. Calculation procedure of the ray-tracing code LHDGauss. This code is incorporated into the AutoAna system [10]. The output file of power deposition profiles will be included as an input file for the TASK3D-a [20]. applied for modelling in this work is at most 1%, which indicates that the magnetic field outside the LCFS is approximated to the vacuum magnetic field. Furthermore, in a case when the neutral beam injection heating is applied, the beam component can modify the pressure profile and the best fit pressure profile can be inconsistent with that expected from kinetic thermal beta. The difference of the refraction of rays in the cases with and without an accounting periphery plasma outside the LCFS should be expected especially in higher density plasmas, although the effect cannot be spatially noticeable under the plasma parameters discussed in this paper due to the current resolution of the ECE diagnostics.
The injection settings for ECH are determined from the ECH setting logs produced each shot. The logs include input power, injection timing, target positions, polarisation states, angles of final steering antennas and Gaussian beam parameters for each EC beam line.
Second, the mode purity is calculated with the 1D fullwave equation along the propagating direction with the cold plasma dielectric tensor, which is summarised in appendix B. Because the contents of each mode excited at the plasma core depend not only on the injection angle and the polarisation states but also on the density gradient and the magnetic shear in the interface region [23,24], the electromagnetic waves with all three components are solved along the propagating direction from the antenna to the target point. In this calculation, the extrapolated equilibrium map is used and all three components of the magnetic field are taken into account to include the density gradient and the magnetic shear. The electric field components perpendicular to the propagating direction are decomposed into two orthogonal components with real and imaginary parts, respectively, in order to express the polarisation state of the electric field. The local dispersion relation using the local n e and the local magnetic field gives the relation between the O-mode and the X-mode electric fields. Each local mode content is determined by checking the orthogonality between the polarisation state of the local electric field and that of the O-or X-mode electric fields, respectively [25], i.e. the local electric field can be decomposed into the orthogonal O-and X-mode electric fields. Each mode content converges to a certain value as the EC wave propagates through the plasma interface region. The excited modes with the fixed O/X mode ratio propagate to the plasma core region, where the wave numbers for both modes are well separated to , where k O and k X are the wave numbers for the O mode and the X mode, respectively, and λ s is the scale length of the magnetic shear.
Third, ray-tracing calculations are executed to obtain the power deposition profiles on each magnetic flux surface labelled by r eff . The injection Gaussian beam is modelled as multiple rays equally distributed with the Gaussian-weighted power fraction for each beam segment. The propagation of each ray is based on the model of geometrical optics with the cold plasma dispersion relation to save the computer resources. The absorption of each ray is calculated according to the quasi-perpendicular weakly relativistic model [26]. The neglected Doppler shift may affect the absorption in the case of oblique propagation of the beams. A possible maximum shift of the absorption location with respect to the magnetic flux surfaces is estimated to be ( for high temperature plasmas. Interactions between the rays in their propagation and absorption processes are not included as are done in the beam-tracing or the quasi-optical beam-tracing code [27][28][29]. For the ray-tracing calculations, local physical quantities, such as ( ) T r e eff and ( ) n r e eff , are obtained by interpolating r eff from the mesh data of the 3D equilibrium mapping, assuming that the T e and n e are constant over a flux surface. The vacuum magnetic field components are also interpolated from the mesh data registered. The absolute values of the power deposition profiles for the O mode and the X mode are obtained by the power ratios between those modes calculated with the 1D full-wave equation described above.
Finally, the ray-tracing calculation results, which are composed of the injection settings, the mode purity, the ray trajectories, the extended equilibrium mapping and the power deposition profiles, are processed as the unified format for the LHD database and registered in the LHD database.
These calculations are automatically executed for each shot and each time slice of the Thomson scattering diagnostic measurement during a pulse, and they are computed in parallel with multiple CPU cores in the AutoAna system. This fast calculation scheme together with comparisons with the Thomson scattering diagnostic measurement or the electron cyclotron emission (ECE) diagnostic enables the seeking of the optimal injection settings for a required power deposition profile on a shot-by-shot basis during ECH experiments in the LHD.
LHDGauss will be introduced into the TASK3D-a [20], the integrated transport analysis suite, as an ECH module so that the calculated power deposition profile is used for the heat transport study for LHD plasmas.

ECH experiments with LHDGauss
In this section, ECE diagnostics signals [30,31] along with the Thomson scattering diagnostics and the absorbed power estimated with the diamagnetic flux measurements [31,32] are compared with the results of LHDGauss for its validation. The dependence of absorbed power on mode purity is also discussed. , , f f f is defined on the plane perpendicular to the radial (R) direction from the centre axis to the outer port #2 centre and the origin of the plane is set at (  figure 2(c), which coincides with the ray-tracing calculation results of on-axis heating. It is observed that the radial widths of the ECE amplitudes are wider than the widths of the absorption profiles by LHDGauss. The widths of the ECE signals are affected significantly by the perpendicular heat transport especially at such a low modulation frequency [33]. The effect will be analysed in the future by using the transport code TASK3D-a, in which LHDGauss will be incorporated as the input for a heat source. Figure 3 shows the absorbed power estimated with the diamagnetic flux measurement and the absorbed power calculated by LHDGauss for each beam line as a function of line-averaged electron density n e,fir measured by the FIR laser interferometer.

Comparison with diamagnetic flux measurement
In these experiments all the injection antennas were targeted toward the vacuum magnetic axis. The absorbed power is experimentally estimated by calculating the change in the time derivative of the plasma's stored energy just before and after MECH on/off timing. The 2-O port 77 GHz injection system, horizontally launched from the outer port #2, has high absorption rates in a relatively high n e,fir region (< × 1.5 10 19 m −3 ) although the rays are refracted and deviated from the magnetic axis for high n e plasmas. The beams from the 2-O port 154 GHz injection systems keep full absorption rates since the rays are not refracted for the range of n e,fir , as shown in figure 3. The one 154 GHz EC beam line (2-OUL) shows the absorption rates above 100%, but these can be attributed to the errors arising from the uncertainty in the estimation of both the absorbed power P abs and the input power P in . Here, '2-OUL' denotes that the end of the waveguide is located on the upper left side of the outer port #2 seen from the outside of the LHD vacuum vessel. It is found that the rays from the 5.5-U port 77 GHz injection system, launched almost vertically from the upper port #5.5, are refracted in high n e plasmas, which causes degradation of the power absorption.
The ray-tracing calculations with LHDGauss show that the absorption rates are almost 100% for the 2-O EC beam lines, which is in rough agreement with the experimental values of the high absorption rates. However, the deviations are large between LHDGauss and the experiment for the 5.5-U EC beam line. One of the reasons is that the effect of multiple reflection of rays is not included in LHDGauss, although the single-pass absorption is weak for the 5.5-U beam line in this experiment. The injection antenna was aligned for on-axis heating in the vacuum magnetic field configuration, which resulted in the refraction of the rays and their deviation from the EC resonance layer. The effect of the multiple reflection can be clearly recognised in figure 4, which shows the modulated components of ECE amplitudes for the 5.5-U and 2-OUR ECH, respectively, as a function of r a / eff 99 . In this discharge, the absorbed power from the 5.5-U ECH is estimated to be almost half that from the 2-OUR ECH, although the two EC beam lines each have almost the same injection power, as high as 1 MW. Again, the injection antennas were aligned for the vacuum magnetic axis. As shown in figure 4, the peaked profile at ∼ r a / 0.55 eff 99 is clearly observed for the 2-OUR ECH, while the almost flat profile is observed for the 5.5-U ECH. This fact suggests that the main part of the injected EC beam from the 5.5-U port does not reach the resonance layer with the single pass, while the injected EC beam from the 2-O port does. This gives rise to a high demand for fine tuning of the 5.5-U ECH injection in order to compensate the refraction effect so as to keep on-axis heating under high n e plasmas. Eventually, since the effect of the multiple reflection is difficult to take into account for ray-tracing calculations due to complicated inner structures of the LHD, the absorption rates tend to be underestimated from the measurement unless single-pass absorption for all rays is achieved.

Dependence of absorbed power on mode purity
The dependence of absorbed power on mode purity was evaluated with LHDGauss, which includes the mode purity  (the outside of the LCFS) shows the extrapolated virtual magnetic flux surfaces. In these calculations, the Cartesian coordinate system is defined by the following unit vectors: , , x y z where k and ϕ denote the unit vectors of the wave number and the toroidal direction, respectively. This coordinate system is independent of the static magnetic field in order to define the polarisation of the injected millimeter waves on the plane perpendicular to k. According to this coordinate system, φ and θ are represented by x . Due to orthogonality of the polarisation, the O/X-mode purity is obtained by are identical and η = σ 0 if they are orthogonal [25]. As shown in figures 5(a1)-(e1), the injected millimetre wave with the given polarisation couples to electromagnetic waves with the initial polarisation remaining in the extremely low-density region of < n 10 e 18 m −3 . The electromagnetic wave has no inherent mode such as the O mode and the X mode, where the two modes are degenerate due to the lowdensity plasma. It is found that the electromagnetic wave gradually couples to the O mode and propagates into the plasma peripheral region. The O mode purity reaches almost unity and it remains inside the LCFS. In the case of nearly perpendicular injection (θ ∼°90 ) from the 5.5-U port, it is expected that the linearly polarised wave (β ∼°0 ) with α same as φ can excite the pure O mode at the peripheral region, where φ ∼°45 . That is why the almost 100% O mode content is obtained by the linearly polarised wave with optimum α ∼°45 set . However, as shown in figures 5(a2)-(e2), the O mode content is obtained at less than unity due to the inappropriate polarisation state of α ∼°25 set . Oscillating electric fields in the propagating direction are observed inside the plasma region due to different wave numbers for each mode.
These calculations may be more useful in a situation in which the magnetic shear and density gradient are larger. Due to the interface region, the exact plasma boundary no longer exists, where the injected millimetre wave couples to the plasma wave. The optimal injection polarisation cannot be determined with the magnetic shear angle at the plasma boundary. The optimal injection polarisation can be finely searched by using this calculation method, which may suggest that elliptical polarisation is needed even for a perpendicular injection case. Figure 6 shows (a) the absorption rate and (b) the mode purity for the 5.5-U ECH as a function of α with β fixed at β° 0 . The absorption rates are experimentally estimated at the ECH turn-on and off timing, which is compared with the ray-tracing calculation results with LHDGauss. The absorption rates calculated with LHDGauss agree well with the experimental data due to inclusion of the mode content analyses. The evaluated absorption rates are represented by a sinusoidal dependence with a suitable maximum value fitted. Small deviations between the mode purity data and the sinusoidal function are considered to be a result of the effect of finite ellipticity necessary for oblique propagation as in the injection from the 5.5-U port.

Optimisation of ECH injection settings
The extension of the high-temperature regime for the LHD was successfully achieved by the utilisation of LHDGauss upgraded with the rapid post-processing magnetic coordinate mapping system. The injection antennas were first aligned for the vacuum magnetic axis. Then, the injection antennas were aligned for the high-temperature plasma, which was based on the ray-tracing calculations of the previous discharges of high-temperature plasmas with the same conditions except for the injection antennas. Figures 7 and 8 show the ray-tracing calculation results (1) before and (2) after tuning the injection antennas of the 5.5-U and the 2-OUR EC beam lines, respectively, for the high-temperature plasma under the experimental n e and T e profiles along with the 3D equilibrium mapping. In figure 7, the rays are projected onto two planes: , where e r denotes the unit vector in the radial direction at the 5.5-U port centre, e tor denotes the unit vector in the toroidal direction at the 5.5-U port centre, and e inj denotes the unit vector in the direction from the injection antenna centre to the target point. The contours of magnetic flux surfaces r a / eff 99 and the electron cyclotron frequencies f ce , the lines of the right-handed wave cutoff frequency f rc , the upper hybrid resonance frequency f uh and the electron cyclotron resonance frequency coincide with 77 GHz, and the first wall is also superimposed on the planes. The resultant power deposition profiles as a function of r a / eff 99 are shown in (c), which indicates that tuning the injection antenna contributes to high power deposition near the magnetic axis, compared to no deposition in the case without the tuning. In figure 8, the rays are projected onto two planes: (a) the plane consisting of two orthogonal basis vectors, are shown in (c), which indicates that tuning the injection antenna contributes to on-axis heating. Figure 9 shows T e and n e profiles as a function of r a / eff 99 in the cases with and without tuning the injection antennas for all the EC beam lines. An almost maximum power as high keV for the case without tuning so far. Besides, a clear increase in T e0 of up to 7.5 keV is observed due to tuning the injection antennas based on the ray-tracing calculations with LHDGauss. This outstanding result is an extension from the conventional high T i plasmas in the LHD [34].

Summary
The ray-tracing code LHDGauss has been upgraded to include the rapid post-processing 3D equilibrium mapping system. By using the interpolated and extrapolated 3D equilibrium  mapping data, the mode contents are determined through solving the 1D full-wave equation along the injection propagating direction from the antenna to the absorption target point. The n e profiles and the magnetic shear affect polarisation of the injected electric fields and determine optimal coupling to plasma waves, thereby changing the absorbed power, which is confirmed in experiments. The optimal coupling can be solved by using this calculation process and this consideration may be more crucial if a magnetic shear and a density gradient are larger, where elliptical polarisation is needed even for a perpendicular beam line.
The ray-tracing calculation results are extensively compared with the LHD experiments. The deposition peak positions calculated with LHDGauss show good agreement with the phase bottoms and the amplitude peaks in the ECE diagnostics in the case of off-axis heating. The on-axis heating cases confirmed with the Thomson scattering diagnostics coincide with the ray-tracing calculation results. The measured absorbed power is comparable to the ray-tracing calculation results if the single-pass absorption rate is almost 100%. The calculated absorption rate tends to be underestimated from the measurement if the single-pass absorption is weak due to lack of the multi-pass absorption effect, which is not implemented in LHDGauss.
The power deposition profiles calculated at each shot and each time slice of the Thomson scattering diagnostics are automatically registered in the LHD database, which contributed to searching for the optimal ECH injection settings during experiments. The adjustment of the injection settings with LHDGauss resulted in a successful increase of T e0 up to 7.5 keV for the high T i plasma, which is an extension of the LHD operational regime. In the near future, higher T e can be expected by optimising not only an injection antenna but also injection polarisation for each EC wave with the help of LHDGauss. 3D equilibrium mapping data The interpolation and extrapolation of the 3D equilibrium mapping data TSMESH that relates ( ( ) ( ) ) ϕ ϕ x r y r z , , , , to r eff are conducted by fitting a 3D bi-quadratic formula, whose coefficients are deduced from minimising the fitting error using the singular value decomposition (SVD) method. Let us suppose that r eff exists in the outside of the LCFS and the extrapolation is executable to the outside of the LCFS although the extrapolation does not have any physical meanings at all. r eff can be expressed in the bi-quadratic formula such that where a 1 = 1, a 2 = x, a 3 = y, a 4 = z, = a x 5 2 , = a y 6 2 , = a z 7 2 , a 8 = xy, a 9 = yz, and a 10 = zx, and ξ j denotes their coefficients, respectively. In a set of sampled mesh points near the object point in the 3D equilibrium mapping space, i.e.   The pseudo-inverse matrix provides a solution to the least squares problem, where the ξ is solved by minimising 2 SVD provides an accurate way to compute the pseudoinverse matrix along with an optimal solution to the least squares problem. SVD is implemented numerically in most computer languages as a function or a module. Interpolation and extrapolation for a given point can be executed in this way. Figure A1 illustrates the original and extrapolated virtual magnetic flux surfaces r a / eff 99 on poloidal planes, for example, at (a) ϕ =°0 , (b) ϕ =°9 , and (c) ϕ =°18 . The r a / eff 99 contours near and outside the LCFS are smoothly generated, although those near the first wall away from the LCFS are not smooth lines.