Current Challenges in the First Principle Quantitative Modelling of the Lower Hybrid Current Drive in

The Lower Hybrid (LH) wave is widely used in existing tokamaks for tailoring current density profile or extending pulse duration to steady-state regimes. Its high efficiency makes it particularly attractive for a fusion reactor, leading to consider it for this purpose in ITER tokamak. Nevertheless, if basics of the LH wave in tokamak plasma are well known, quantitative modeling of experimental observations based on first principles remains a highly challenging exercise, despite considerable numerical efforts achieved so far. In this context, a rigorous methodology must be carried out in the simulations to identify the minimum number of physical mechanisms that must be considered to reproduce experimental shot to shot observations and also scalings (density, power spectrum). Based on recent simulations carried out for EAST, Alcator C-Mod and Tore Supra tokamaks, the state of the art in LH modeling is reviewed. The capability of fast electron bremsstrahlung, internal inductance li and LH driven current at zero loop voltage to constrain all together LH simulations is discussed, as well as the needs of further improvements (diagnostics, codes, LH model), for robust interpretative and predictive simulations.


1.Introduction
With the highest known current drive (CD) efficiency in tokamaks and off-axis power absorption, the rf wave at the Lower Hybrid (LH) frequency is particularly attractive for current profile shaping or steady-state operation [1].In this context, first principle modeling has started long time ago and is still active after more than 30 years of research.Thanks to the challenges concerning ITER physics performances by appropriate current profile tailoring, an active international modeling activity has been carried out during the last five years, involving quantitative comparisons between simulations and experiments for several tokamaks (Tore Supra, Alcator C-Mod, EAST) [2,3].In this paper, a review of main achievements is presented, which highlights the state of the art in this field of physics.While modeling tools are able to reproduce quantitatively observations in low density plasmas, simulation becomes challenging when LH wave absorption becomes weak, i.e. for plasmas at high density and low temperature.In this regime characterized by a rather large spectral gap, the level of agreement between simulations and experimental results may be significantly improved by a broadening of the launched power spectrum at the separatrix.The potential influence of edge plasma physics on core absorption is one the many challenges to be tackled concerning LH current drive in tokamak, in particular for ITER.

Numerical tools used by LH current drive first principle modeling
The aim of first principles modeling is to build a consistent interpretative framework of LH experiments with reliable predictive capabilities for existing tokamaks and ITER reactor.Such an objective represents a considerable challenge because of the multiple physical mechanisms involved concerning propagation and absorption.To this day, there is still a remarkable contrast between the reproducibility of LH experiments and the sensitivity of LH modelling.Multi-machine modeling using well benchmarked numerical tools with quantitative comparison against experimental observations is the reference path to reach this ambitious goal.
As shown in Fig. 1, the general scheme for first principle modeling consists to couple different models describing successively wave propagation from the antenna in the Scrape-Off Layer (SOL), while wave propagation and absorption of the LH wave inside the separatrix are determined by solving numerically Maxwell and Boltzman equations in the geometrical optics and Fokker-Planck approximations respectively [2].The standard tools used for this purpose are C3PO/LUKE or GENRAYCQL3D raytracing/Fokker-Planck codes [4,5,6,7].In presence of a tail of fast electrons that are pulled out from the thermal bulk by resonant Landau interaction with the LH wave to high kinetic energies, nonthermal bremsstrahlung in the hard Xray (HXR) photon energy range provides a natural diagnostic with an almost direct insight on the electron momentum distribution function built-up at different radial positions, an invaluable diagnostic to validate models and numerical tools.The R5-X2 code coupled to LUKE calculates the fast electron bremsstrahlung (FEB) and an analogous calculation is performed inside CQL3D itself, using the same numerical kernel for the quantum-relativistic cross-sections [8].
All the concepts for first principle modeling LH wave current drive were already established about 30 years ago.Since, progresses have concerned principally more realistic descriptions of the LH wave and electron dynamics in momentum and radial spaces, allowing systematic quantitative comparisons between simulations and experiments, which represents itself a significant achievement.In the meantime, powerful toroidal full-wave codes have been developed (TORLH, LHEAF) [9,10,11].Although both are the same full-wave code for LH waves, the two codes took different sort of complimentary approaches.TORLH started from TORIC, which uses a spectral basis and solved a dense matrix problem, while LHEAF uses the finite element method and solves a sparse matrix problem.Both codes have advantages and issues.For the LH wave, the electron Landau damping is the only term which requires a dense matrix and even TORLH simplifies a dispersion relation to eliminate hot plasma branch.This fact has been exploited in LHEAF, making the code becomes more efficient in terms of computer resource.However, the approach in LHEAF eliminates some of physics such as a mode conversion at the LH resonance, which could be of interest to investigate in some specific situations.On the other hand, the spectral approach does not allow a realistic description of the SOL while some further improvements must be considered for the launcher implementation in TORLH.It is important to recall that TORLH has been used to validate a posteriori some critical assumptions considered in ray-tracing calculations, which is particularly important for present day simulations [9].A systematic use of full-wave codes for first principle modeling is likely one the challenges in this domain, thanks to the very large computational effort, whatever the considered method.
The first step in LH current drive modeling is to calculate the power spectrum that must be considered at the plasma separatrix, by performing of FFT of the wave electric pattern along the toroidal direction (z).This is done by wave coupling codes like ALOHA [12], considering the antenna structure as arrays of square waveguides along the toroidal and poloidal directions.Depending upon the waveguide phasing along the toroidal direction, the LH power spectrum is made of a set of lobes characterized by a Gaussian-like shape, each of them centered at a toroidal refraction number N z , with a width ∆N z resulting from a FFT of the antenna size in the toroidal direction.Uniform variations of density in front of the antenna by any process (ponderomotive effects at high input power, gas feed,…) do not modify spectral lobe positions, but their relative heights, changing the overall directivity, i.e. the relative fraction of power that generates co-or counter-current in the plasma [13].The determination of the directivity has a large impact in code predictions of the total LH current, especially for antennas who turns out to be very sensitive to edge density variations like at 2.45 GHz in EAST tokamak.In this case, the consistency between edge density deduced from reflection coefficient (RC) and measured directly by probes in front of the antenna is essential to guarantee the robustness of the launched power spectrum.For the antenna at 2.45 GHz in EAST tokamak, a large discrepancy between RC and density measurement is observed as shown in Fig. 2, a problem already encountered in Tore Supra LH experiments [14].
If the SOL is turbulent, depending upon the spatial characteristics of the density fluctuations, lobe positions may change randomly, as well as their amplitude [15].This approach is providing a a posteriori interpretation of the tail LH model, which was introduced to explain successfully HXR profiles in Tore Supra [16].The transfer of power from the main lobe to the tail is a free parameter, as illustrated in Fig. 3, and significant effects on the predicted current density and HXR profiles require at least that P tail /P tot > 0.5.Beyond this threshold, predicted results are weakly sensitive to this ratio, making the model rather robust.Ray-tracing calculations are used for describing LH wave propagation that is separated from the absorption process, as weak damping approximation is always satisfied.The well collimated beam and WKB approximations are applicable to a large domain in the plasma, except at cutoff and caustics, which may represent an issue when ray length is long in weak absorption regime.However, thanks to full-wave calculations, errors are likely benign for cut-off near the plasma edge, while, at caustics, the impact of this problem on overall simulations remains more open [17].Since ray-tracing captures toroidal refraction from B-field inhomogeneity, a natural mechanism in tokamaks due to toroidal curvature, it is an appropriate tool for bridging the spectral gap leading to wave absorption, while for N ||0 corresponding to the main lobe, the plasma is in principle transparent.However, this mechanism is not necessarily the single one that could contribute to LH wave absorption, and the difficulty in first principle modeling is to identify the possible existence of additional external physical mechanisms for bridging the spectral gap.
The condition of absorption corresponds to large N || upshift combined with inward propagation, which is difficult to satisfy for plasmas with circular poloidal cross-section or large aspect ratio tokamaks, while small aspect ratio diverted plasmas are much more favorable.This is an important difference between magnetic configurations.Indeed, if ray stochasticity develops before absorption, LH simulations are becoming very sensitive to any small variations of the toroidal MHD equilibrium or initial ray condition, a problem never observed experimentally [16].

(Top) Waveform of the reflection coefficient (RC) and the measured density in front of the 2.45 and 4.6 GHz antennas of the EAST tokamak; (down) RC variation with density in front of the antenna and corresponding LH power spectra (insert) as calculated by the wave coupling code
ALOHA [12].For the 2.45 GHz antenna, the measured RC level is corresponding to a density of ~4x10 18 m -3 while direct density measurement is 4 times lower.Furthermore, it is shown that the power spectrum of the 2.45 GHz antenna is significantly modified by changing the density (red curve in the insert: ~4x10 18 m -3 , blue and black curves: ~7x10 17 -1x10 18 m -3 ).For the 4.6 GHz antenna, RC level and density measurement are consistent, leading to a robust estimate of the launched power spectrum.
This suggests that additional mechanisms may come into play to complement weak toroidal refraction, and bridge the spectral gap before wave penetration in the plasma core.To tackle this difficulty, a large number of rays with a small numerical width is considered using GENRAY/CQL3D calculations, assuming that numerical solution converges towards the physical one as the number of rays is large [18].For C3PO/LUKE, each lobe in the power spectrum is considered locally as a plane wave, even for the tail LH model, with a spectral width given by the Fourier transform of the initial wave electric field pattern.With this approach, all rays are well separated from each other in phase space, a condition for a correct calculation the LH diffusion coefficient by simply adding all ray contributions [2,16].
The LH wave absorption is calculated by a 3-D bounce-averaged relativistic Fokker-Planck solver using a linearized collision operator (either in LUKE or CQL3D) that determines the distortion of the electron distribution functions under the influence of the LH wave electric field.The corresponding rf diffusion operator describing the resonant Landau absorption (RLA) is determined in LUKE and CQL3D codes from the relativistic bounce-averaged Kennel-Engelman-Lerche formula for plane waves, taking into account of the exact LH wave polarization (slow and fast modes).Concerning the power decrement for each ray, both CQL3D and LUKE have similar approaches, taking into account of multiple passes of the same ray through an incremental plasma volume between neighboring magnetic flux surfaces.The self-consistency between the electron distribution and the rf diffusion is achieved until all LH power is transferred to electrons for each ray.The quasilinear calculation is performed while other processes may contribute, like parasitic non-resonant collisional absorption (NRCA) which may be important in very cold regions of the plasma, and anomalous fast electron transport, an ad-hoc free parameter in simulations [19].
It is important to recall that RLA calculations are valid for circulating and trapped electrons far from low order rational magnetic flux surfaces (not valid for CD in islands).The calculation of the power transferred from the waves to the plasma is based on the narrow beam approximation in the poloidal direction, i.e. plasma must be almost homogeneous transverse to the beam, a condition which is sometimes not well fulfilled.Besides, ray-tracing does not provide any prescription for ∆N || evolution in the plasma.But it turns out that results are weakly sensitive to ∆N || along the ray trajectory as long as ∆N || << N || corresponding to the well collimated beam approximation [6].
Finally, the ultimate step in first principle modelling consists in determining moments of the electron distribution function, the so-called synthetic diagnostics.They allow direct comparison between simulations and experiments, and are crucial steps in the simulation process.The convolution of the physical signal with the detector's response leads to an additional complexity in the physical interpretation.So, it is important to recall that HXR  tomography in a poloidal cross-section of the plasma is indicative of the LH power absorption profile instead of LH driven current density profile [14].Furthermore, the standard effective charge Z eff cannot be used readily for HXR calculations, except if the plasma contains only fully stripped low Z impurities.For coming experiments with tungsten walls, the use of impurity transport codes for HXR calculations must be systematic [8].Besides HXR measurements, internal inductance, which is a normalized integral moment of the current density profile, is supposed to give also valuable indications on the current density profile [3].However, combined with the plasma current, it cannot constrain unambiguously simulations as illustrated in Fig. 4 for a set of EAST tokamak discharges.Only quantitative measurements that combine simultaneously dynamics in momentum and configuration spaces like HXR are more selective between models [2].Time modulation of the LH power or the antenna phasing may represent also a challenging constraint for first principle modeling.

A selection of results
Quantitative "snapshot" simulations of a single shot in almost steady-state conditions thanks to the resistive current diffusion have been performed several times for Tore Supra, Alcator C-Mod and EAST tokamaks.For low density discharges, a good quantitative agreement is found between simulations and observations concerning HXR line-integrated profiles and the plasma current, as shown for Tore Supra in Fig. 5 [14].This represents an important validation of the whole chain of calculations.For plasmas of higher densities, large discrepancies are observed between FEB profiles and simulations, using the standard LH model (power spectrum with a quiescent SOL).In order to recover a good agreement, the introduction of a fast fluctuating power spectrum (tail LH model) with respect to the fast electron slowing-down time must be introduced in LH-driven Tore Supra plasmas [16].This is similar for Alcator C-Mod as shown in Fig. 6, the tail LH model allowing to recover simultaneously the bell-shape HXR line-integrated profiles observed up to n e ≤ 10 +20 m -3 and the total plasma current [2].Attempt to reproduce the decrease of the internal inductance as the density rises at V loop = 0V is a challenging exercise.This study has been carried out for several sets of EAST tokamak LH-driven discharges, using interfero-polarimetry to enhance quality of the toroidal MHD equilibrium reconstruction by EFIT code [3].As shown in Fig. 7, the decrease of l i is well predicted using C3PO/LUKE using the standard LH model, but the variation of the internal inductance is too large as compared to observations, indicating that the calculated current profile is likely too off-axis.The relative agreement between C3PO/LUKE and GENRAY/CQL3D predictions is remarkable with the standard LH model, except that the current predicted by GENRAY/CQL3D is higher than measurements, while for C3PO/LUKE, its absolute level is very close to the experimental value in the lowest density case.The tail LH model gives conversely the correct current level at all densities, but profiles are too peaked, leading to l i values higher than those deduced from magnetic measurements.From calculations, it could be deduced that LH power absorption is less off-axis than the standard LH model, and more off-axis for the tail one, which highlights the importance of the LH power spectrum considered for ray-tracing calculations in prediction of the current density profile.Since l i is very sensitive to small variations of the current density at mid-radius and outside, the use of diagnostics that provide additional constraints in this region of the plasma is essential.The multi-code analysis illustrates the challenge of current density calculation from first principles, and the importance of the launched power spectrum at the separatrix on general conclusions.Its direct experimental determination will be likely one of the most important challenge to decide on the relative influence between core and edge LH physics on the overall simulations.This problem is crucial for ITER LH current drive predictions, as already pointed out [2].
Scaling of the current drive efficiency with electron density has been found to be highly anomalous in Alcator C-Mod, considering the HXR intensity of central chords as a proxy [20].From GENRAY/CQL3D calculations, the sharp decrease above 10 +20 m -3 is ascribed to parasitic absorption of the LH wave in the cold and dense SOL.The LH power available for current drive in the core plasma is then dramatically reduced.Simulations of HXR emission at 1.3x10 +20 m -3 with C3PO/LUKE/R5-X2 codes have shown a good agreement with observations, without invoking NRCA at the plasma edge as show in Fig. 8 [2].
For this calculation, only 5% of the total current at this density is driven by the LH wave, which leads to questioning the non-inductive nature of the discharge [2].For densities lower than 10 +20 m -3 , as foreseen in ITER, the experimental HXR scaling with n e is grossly reproduced with the standard or tail LH model, either by C3PO/LUKE and GENRAY/CQL3D codes, as shown in Fig. 8.This scaling is also well reproduced in Tore Supra up to 5.5x10 +19 m -3 as shown in Fig. 9, but the quality of the agreement between modelling and experiment is greatly improved using the tail LH model [16].

Conclusions
For small spectral gap regimes, the existing first principle modeling tools are able to describe observations quantitatively, which represents a remarkable achievement.However, large spectral gap regimes remain still challenging to model, and it is difficult to decide if the discrepancy arises from a failure in the chain of modeling or from physical processes that could be missed in LH current drive modeling at the edge or in the core of the plasma.In this context, one of the most important challenge is to determine experimentally the characteristics of the power spectrum as the LH propagates through the SOL, considering the large impact for LH current drive in ITER.But, numerous other investigations must be carried out to describe simultaneously HXR profiles and l i at V loop = 0V and to improve the consistency between measured density in front of the antenna and the RC level.Besides, simulations must be performed on a much larger database, to estimate the robustness of LH simulations.These short terms objectives must not hide the numerous studies that must be revisited, like Ohmic/LH synergy, each of them being able to put larger constraints on modeling, in particular when quasilinear effects become very important.Finally, on a long term, the fullwave calculations must likely become a standard, in replacement of ray-tracing calculations, in such a way all intrinsic limitations may be removed and new synthetic diagnostics must be developed.

Figure 1 .
Figure 1.General scheme for first principle modelling for rf heating and current drive in tokamaks.

Figure 2 .
Figure 2. (Top) Waveform of the reflection coefficient (RC) and the measured density in front of the 2.45 and 4.6 GHz antennas of the EAST tokamak; (down) RC variation with density in front of the antenna and corresponding LH power spectra (insert) as calculated by the wave coupling codeALOHA[12].For the 2.45 GHz antenna, the measured RC level is corresponding to a density of ~4x10 18 m -3 while direct density measurement is 4 times lower.Furthermore, it is shown that the power spectrum of the 2.45 GHz antenna is significantly modified by changing the density (red curve in the insert: ~4x10 18 m -3 , blue and black curves: ~7x1017 -1x10 18 m -3 ).For the 4.6 GHz antenna, RC level and density measurement are consistent, leading to a robust estimate of the launched power spectrum.

Figure 3 .
Figure 3. Standard and tail power spectra for a main lobe at N ||0 .The tail spectrum is made of 5 secondary lobes up to N ||max , which is generally chosen as 6.5/√T e0[16].

Figure 4 .
Figure 4. (Top) waveforms of LH-driven EAST discharges at V loop = 0 V, where the internal inductance l i is ranging from 1.0 (green), 1.1 (blue), 1.2 (red) to 1.4 (black); (Down) (a) and (b) are peaked and off-axis current density profiles which reproduce both the range of observed l i : 1.0 (black), 1.1 (green), 1.2 (blue) to 1.4 (red).Curves in (c) and (d) are corresponding normalized lineintegrated HXR profiles assuming proportionality to the current density, while dots are experimental measurements.Because of their statistical noises, both hollow (a) and peaked (b) current density profiles are consistent the normalized HXR measurements, whatever the l i value.

Figure 8 .Figure 9 .
Figure 8. Scaling of the central chord HXR intensity as function of density for Alcator C-Mod.Full line: GENRAY-CQL3D (standard LH model without NRCA) dashed line: GENRAY-CQL3D (standard LH model with NRCA), yellow points: C3PO/LUKE (standard LH model without NRCA), red points: C3PO/LUKE (tail LH model without NRCA).Other points are experiments[20,2] Figure9.Scaling of the central chord HXR intensity as function of density for Tore Supra, using C3PO/LUKE with tail LH model without NRCA.With the standard LH model, the tendency is similar, but predictions are noisy, because of the ray stochasticity[16].