Plasma and Fusion Research: Regular Articles

Volume 11, 2403098 (2016)

Comparison between Full Wave and Ray-Tracing Calculations toExamine Scenarios for Electron Bernstein Wave Heating in LHD∗)
Hiroe IGAMI, Atsushi FUKUYAMA1) , Hiroshi IDEI2) , Kazunobu NAGASAKI3) , Yuki GOTO4) ,Shin KUBO, Takashi SHIMOZUMA, Yasuo YOSHIMURA, Hiromi TAKAHASHI,Toru I. TSUJIMURA and Ryohei MAKINONational Institute for Fusion Science, 322-6 Oroshi-cho Toki 509-5292, Japan
Department of Nuclear Engineering, Kyoto University, Nishikyo-ku, Kyoto 615-8530, Japan2)
Research Institute for Applied Mechanics, Kyushu University, 6-1 Kasuga-koen, Kasuga 816-8580, Japan3)
Institute of Advanced Energy, Kyoto University, Uji 611-0011, Japan4)
Department of Energy Engineering and Science, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8603, Japan1)(Received 30 November 2015 / Accepted 18 April 2016)
With the use of the density scale length and the magnetic field strength around the O-X mode conversion
region in the LHD for 77 GHz EC waves, the full wave calculations by the TASK/WF2D code have been performed. The wave patterns obtained by full wave calculation have been compared with the trajectories obtainedby multi ray tracing calculations. For the X-mode and the O-mode reflected near each cutoﬀ, wave patterns obtained by full wave calculations and trajectories obtained by multi ray-tracing correspond to each other well. On
the contrary, trajectories of the rays restarted at the high field side of the evanescent region between the plasma
and the left handed cutoﬀs do not cover fully the wave patterns after the O-X mode conversion obtained by fullwave calculations. Improvement of the manner to restart the ray-tracing is required. This improvement can provide a more precise view of the propagation of the EBWs and may contribute to examining scenarios for electron
Bernstein wave heating.
c 2016 The Japan Society of Plasma Science and Nuclear Fusion Research

Keywords: electron cyclotron resonance heating, electron Bernstein wave, full wave calculation, ray-tracing
DOI: 10.1585/pfr.11.2403098

1. Introduction
Electron cyclotron resonance heating and current
drive (ECRH and ECCD) in over dense plasmas wherethe electron density is more than the cutoﬀ density of theelectron cyclotron (EC) wave are available with use of theelectron Bernstein waves (EBWs). Since the EBWs are
the electrostatic waves, they should be excited via the linear mode conversion process inside the plasma by launching the electromagnetic (EM) waves from the vacuum region [1]. To excite the EBWs eﬃciently when the density scale length in the “mode conversion region” is longenough compared to the wavelength, as in the over dense
plasma generated in the Large helical device (LHD), theEM waves should be injected with an appropriate polarization to couple with the O-mode at the plasma boundary
toward an appropriate direction so that the parallel refractive index N// consists with [2] or is very close to the op√
timum value N//opt = Ωce /(ω + Ωce ) at the plasma cutoﬀ(PC), where Ωce is the cyclotron angular frequency and ω is
the wave angular frequency. The injected O-mode couples
with the X-mode that propagates toward the upper hybrid
resonance (UHR) and finally couples with the EBW.
author’s e-mail: igami@LHD.nifs.ac.jp∗) This article is based on the presentation at the 25th International Toki
Conference (ITC25).

In the LHD, ECRH in over dense plasmas by EBWs
excited via the O-X-EBW mode conversion process mentioned above was performed [3]. In the analyses of theEBW heating experiment, the multi ray-tracing taking into
account the finite beam width is used. However, the raytracing cannot provide a precise view of the wave propagation around the mode conversion region because the WKB
approximation is not valid there. When N// is not strictly
equal to N//opt at the PC the ray-tracing shows that the injected O-mode is just reflected in front of the evanescentregion even if a part of the power of the injected O-mode
can tunnel through the evanescent region between the PC
and the left handed cutoﬀ (LC) and can couple with the Xmode with the eﬃciency as the analytic form written below
suggests [4].
T OX = exp{−πk0 Ln (β/2)1/2
× [2(1 + β)(N// − N//opt )2 + Nv2 ]},(1)

where k0 = ω/c, c is the velocity of light in vacuum, Ln isthe scale length of the density gradient, Nv is the refractive
index perpendicular to the magnetic field and the density
gradient.
To investigate the wave propagation and absorption after the tunneling through the evanescent region, an “expe-2403098-1c 2016 The Japan Society of Plasma
Science and Nuclear Fusion Research

Plasma and Fusion Research: Regular Articles

Volume 11, 2403098 (2016)

dient manner” to restart the ray-tracing calculation at theedge of the higher density side of the evanescent region
has been used to analyze the wave propagation in the overdense region [5]. However, to describe the wave propagation around the mode conversion region precisely thefull wave analysis that solves the Maxwell’s equation directly is required with taking into account the non-uniform
plasma profiles and the distribution of the wave electricfield. We are developing a numerical code including theTASK/WF2D code that solves the Maxwell’s equation by
finite elements method (FEM) in two-dimensional space[6]. This numerical code can describe the wave propagation across the evanescent region with taking into accountthe non-uniform plasma profiles and arbitrary distributionof the electric field of the incident wave. In this paper,comparisons between the full wave calculations obtainedby the TASK/WF2D code around the evanescent region
and the ray-tracing calculations with use of the “expedient manner” to check its validity is reported. In section 2,the systems for the full wave calculations (TASK/WF2D)
and the ray-tracing calculations are explained. In section
3, comparison of the calculation results by full wave andray-tracing calculations with use of the parameters aroundthe O-X mode conversion region in LHD are shown. Insection 4, we summarize the contents.

2. Systems for the Full Wave and RayTracing Calculations
TASK/WF2D solves the wave equation as the boundary problem∇ × ∇ × E − k20

ε · E = iωμ0 jext ,(2)value of the electric field is given at x = xw . In the current
version of TASK/WF2D, plasma parameters that change
linearly along the z and x directions can be given. The
permittivity tensor of the collisional cold plasma is used.
The electric field as the boundary value is given basedon the electric field at the aperture plane of the rectangular waveguide whose inner width is a. At the aperture,the electric field strength is given as the basic rectangular waveguide mode E = E0 sin(πl/a), where l is the distance from (xw , z1 ) along the z direction shown in Fig. 1and a is the length of the inner side of the waveguide.
When the injection angle between the wave vector and thez direction is θ, the electric field strength along the z direction at x = xw is given with taking into account thephase delay of the plane wave front so that the phases φ1at (xw , z1 ) and φ2 at (xw , z2 ) should be set as input parameters so that dφ = φ2 − φ1 = k0 (z2 − z1 ) cos θ. In the region z1 ≤ z ≤ z2 at x = xw the y and the z componentsof the electric field are given as Ey = E(sin ψ − ie cos ψ),and Ez = E(cos ψ + ie sin ψ) with use of real value input parameters E, e, and ψ to determine relative electricfield strength and phase. To show the wave propagation
across the evanescent region clearly without the reflection
at the conductive wall, the wave absorber whose thickness is wdamp given as an input parameter can be put inalong the calculation boundary without the region between(xw , zdmin ) and (xw , zdmax ). Here zdmin and zdmax are given
as input parameters.
In the ray-tracing calculations, the same magnetic
field and the electron density profile as used in the fullwave calculations by the TASK/WF2D code are adopted.
The dispersion equation of the cold collisionless plasmaaN⊥4 − bN⊥2 + c = 0

where ε is the dielectric tensor, E is the electric field vector, μ0 is the vacuum permeability, jext is the external current vector. Figure 1 shows a schematic view of the system of TASK/WF2D used in the analysis here. The calculation area is defined with input parameters xc , xw , zmin ,and zmax . The electric field is calculated at each meshed
position along the x and the z directions. The boundary

is used, where a = S , b = RL + PS − (P + S )N//2 ,c = P(R−N//2 )(L−N//2 ), S = 1−α/(1+β2 ), R = 1−α/(1−β),
L = 1 − α/(1 + β), P = 1 − α, α = (ωp /ω)2 , β = Ωce /ω,ωp is the plasma angular frequency. As mentioned in section 1, an “expedient manner” is used to restart the ray-Fig. 1 Schematic view of calculation system for TASK/WF2D.

Fig. 2 Schematic view of the “expedient manner” to restart theray-tracing at the high field side across the evanescent region.

2403098-2(3)

Plasma and Fusion Research: Regular Articles

Volume 11, 2403098 (2016)

tracing after the O-mode tunnels through the evanescentregion and couples with the X-mode in the high density
side. Figure 2 shows the schematic view of the “expedient manner”. When the injected O-mode is reflected near
the evanescent region between the PC and the LC, we determine the reference point where the perpendicular component of the refractive index N⊥ becomes the minimum.
The wave vector at the reference point is taken to be a conserved value ks and look for the point where ks is the solution of eq. (3) at the high density side across the evanescentregion along the direction of the density gradient.

3. Comparison of the Calculation Results by Full Wave and Ray-Tracing
Calculations
In the calculations, the wave frequency 77 GHz that isused for ECRH in the LHD is selected. Since the typicaldensity scale length normalized by the vacuum wavelengthof the 77 GHz EC wave Ln /λ0 = 53 and β = Ωce /ω =0.615 at the plasma cutoﬀ in discharges when the ECRHby EBWs excited via the O-X-EBW mode conversion isperformed, we adopt a simple linear electron density profile changes only along the x direction with Ln /λ0 = 53
and a uniform magnetic field to be β = 0.615 parallel tothe z direction. In TASK/WF2D calculation, we assume
the case where the length of the inner side of the waveguide, a = 0.1 m. The selected mesh spacing of the fullwave calculation is 4 mm.
Note that in the real LHD configuration, the electron
density and magnetic field change three-dimensionally andthe EM waves are launched as a Gaussian beam. However,
it is worth to examining the propagation of the waves thathave a finite beam width to verify the validity of the “expedient manner” to restart the ray-tracing.
Figure 3 shows a wave pattern as a color map of thereal x component of the electric field strength and a distribution map of the absorbed power by collisional damping
as a gray-scaled image near the O-X mode conversion region when the O-mode dominant EM waves are injectedwith the injection angle θinj = θopt = arccos(N//opt ) fromlower right side. A part of the injected EM waves reaches
the high density side of the PC and propagates backward
toward the UHR. The distribution of the absorbed powersuggests this backward wave is fully absorbed by collisional damping near the UHR. Although the launching angle θinj , is set to be θopt , split wave patterns that might bereflected as the O-mode near the PC are shown. Since theradiation pattern of the waves launched from the waveguide aperture is not strictly the same as that of the planewaves, the broadening of the wave pattern might cause deviation of N// from N//opt and reflection near the PC.
With taking into account the broadening of thelaunched wave we performed multi ray-tracing calculations with θinj = θopt at the beam center, θinj = θopt ± 1.5
degrees at each edge of the beam, and plotted the trajec-Fig. 3 Color map of the real x component of the electric fieldstrength (left) and gray-scaled image of the absorbed
power by collisional damping (right) when θinj = θopt .

Fig. 4 (From the left) color maps of the real x, y, z componentsof the electric field strength and gray-scaled image of theabsorbed power when θinj = θopt . Multi ray trajectories of
the reflected X-mode, reflected O-mode, and mode converted X-mode are plotted over the maps. At beam centerθinj = θopt , and at each edge of the beam θinj = θopt ± 1.5.

tories over the wave patterns in Fig. 4, where the real x,y, z components of the electric field strength are shown ascolor maps. Trajectories of the O-mode reflected near thePC are plotted except the central ray of N// = N//opt . Trajectories of the reflected O-mode whose N// deviates from
N//opt a little, well correspond to the split wave pattern of
the reflected O-mode mentioned above. Wave patterns thatmight correspond to the X-mode reflected near the RC are
also shown. The multi-ray trajectories of the reflected Xmode correspond well to the wave patterns. On the otherhand the trajectories of the rays restarted with the “expedient manner” explained in Fig. 2 do not fully cover thewave patterns of the backward waves after the O-X mode
conversion that is obtained from TASK/WF2D calculation.
For the cases of N// = N//opt + 0.02 and N// = N//opt +0.1, full wave and multi ray tracing calculations were per-2403098-3

Plasma and Fusion Research: Regular Articles

Volume 11, 2403098 (2016)

Fig. 5 Similar plots to Fig. 4 when θinj = arccos(N//opt + 0.02).

Fig. 7 Schematic view of the improved idea of the “expedientmanner” to restart the ray-tracing.

to take the appropriate domain of the asymptotic line of thetrajectory drawn in Fig. 7.

4. SummaryFig. 6 Similar plots to Fig. 4. when θinj = arccos(N//opt + 0.1).

formed in the same manner. As shown in Fig. 5, the deviation of N// from N//opt by 0.02 allows eﬃcient O-X mode
conversion while the deviation of N// from N//opt by 0.1 allows a very poor O-X mode conversion rate as shown inFig. 6. Trajectories of the ray of the reflected waves withtaking into account the beam broadening also correspondwell to the wave pattern obtained by full wave calculationsin these cases. Since the N// at the beam center deviatesfrom N//opt the reflection of the O-mode occurs at the beam
center and the wave patterns are not split as in the case ofN// = N//opt . Trajectories of the rays after the O-X mode
conversion do not cover the wave patterns obtained by fullwave calculations fully also.
As explained in Fig. 2, the restart point of the raytracing is determined as a point shifted along the densitygradient (the x direction in this case) from the reference
point. However, the full wave calculations suggest the tunneling waves actually proceed obliquely upward. It might
be better to place the restart point with taking into accountthe propagation toward the direction perpendicular to thedirection of the density gradient (the z direction in thiscase) as shown in Fig. 7. However, verifications with thefull wave and ray tracing calculations are required for how

For the analysis of the ECRH by EBWs excited via
the O-X-EBW mode conversion process performed in theLHD, a numerical code including the TASK/WF2D code isunder development. As the first step, full wave calculations
by TASK/WF2D with the use of the density scale lengthand the magnetic field strength around the O-X mode conversion region in the LHD for 77 GHz EC waves have been
performed and compared with the multi ray tracing calculations. For the O-mode and X-mode reflected near theircutoﬀs, the wave patterns obtained by full wave calculations and the trajectories obtained by multi ray-tracing calcuations correspond to each other well if the broadeningof the injected beam is taken into account. On the other
hand trajectories of the multi rays after the O-X mode conversion do not cover the wave patterns fully. It has beensuggested that it might be better to place the restart pointwith taking into account the propagation toward the direction perpendicular to the direction of the density gradient.
Verifications with the full wave and ray tracing calculations will be required for the improvement of the “expedient manner” to determine the restart point. This improvement can provide a more precise view of the propagationof the EBWs and may contribute to examining scenarios
for electron Bernstein wave heating.
[1] H.P. Laqua, Plasma Phys. Control. Fusion 49, R1 (2007).
[2] J. Preinhaelter and V. Kopecky, J. Plasma Phys. 10, 1(1973).
[3] H. Igami et al., Plasma Fusion Res. 7, 2402110 (2012).
[4] E. Mjølhus, J. Plasma Phys. 31, 7 (1984).
[5] H. Igami et al., Plasma Sci. Technol. 11, 430 (2009).
[6] A. Fukuyama, H. Igami and H. Idei, 41st European Physical Society Conference on Plasma Physics (EPS 2014)
P4.014 (2014).

2403098-4