Plasma and Fusion Research: Regular Articles

Volume 11, 2403070 (2016)

Kinetic Full Wave Analysis of Electron Cyclotron Wave ModeConversion in Tokamak Plasmas∗)
Shabbir A. KHAN1,2) , Atsushi FUKUYAMA1) , Hiroe IGAMI3) and Hiroshi IDEI4)
1)

Department of Nuclear Engineering, Kyoto University, Kyoto 615-8540, Japan
2)
National Center for Physics, QAU Campus, Islamabad 45320, Pakistan
3)
National Institute for Fusion Science, Toki 509-5292, Japan
4)
Research Institute for Applied Mathematics, Kyushu University, Kasuga, Fukuoka 816-8580, Japan(Received 2 December 2015 / Accepted 18 March 2016)
By employing the integral formulation of dielectric tensor, we have introduced the kinetic full wave analysisof ordinary-extraordinary-Bernstein (O-X-B) mode conversion in tokamak plasma in one dimension on the basis
of TASK/W1 code using the finite element method. The boundary value problem of Maxwell’s equation is solved
and the finite Larmor radius eﬀects are represented by integral formulation. The O-X-B mode conversion and theabsorption near the cyclotron harmonic resonance are successfully described.
c 2016 The Japan Society of Plasma Science and Nuclear Fusion Research

Keywords: dielectric tensor, full wave analysis, O-X-B mode conversion, Larmor radius, Bernstein wave
DOI: 10.1585/pfr.11.2403070

1. Introduction
The electromagnetic waves in the electron cyclotron(EC) range of frequencies have been successfully applied
to the electron cyclotron heating (ECH) and current drive(ECCD) in fusion devices. EC waves are extremely useful due to the fact that they can be launched far from theplasma and have features of highly localized and controllable heating capabilities.
The remarkable advantages of ECH and ECCD inhigh density plasma core have led great interest in theoretical and experimental investigations of EC waves in pastdecades [1–5]. Electromagnetic (EM) plasma waves suchas ordinary (O) or extraordinary (X) waves suﬀer the wave
cutoﬀs in high-density regions. This is particularly noticed
in spherical tokamaks (STs) operated in high-β regimes
generally, in which the usual O- and X-modes are cut oﬀ.
In this case, electron Bernstein waves (EBWs) seem tobe the option which provide the penetration to the highdensity region.
Bernstein waves need to be excited externally throughmode conversion of injected EM waves owing to hightemperature in fusion plasmas. Well known schemes ofmode conversion include high field side launch of wave,
direct eXtraordinary Bernstein (X-B) mode conversion and
Ordinary-eXtraordinary-Bernstein (O-X-B) mode conversion. Among them, the O-X-B mode conversion technique [6] is the most promising of ECH/ECCD in overdense plasmas both in tokamaks and stellarators [7–9]. At
specific launching angle, an O-mode wave converts to Xmode at a cutoﬀ frequency. After this conversion, the wave
author’s e-mail: shabbir.khan.36m@st.kyoto-u.ac.jp∗) This article is based on the presentation at the 25th International Toki
Conference (ITC25).

propagate backwards to the plasma edge till it encounters
the upper hybrid resonance (UHR) layer. Since electron
Bernstein wave is a quasi-electrostatic mode required to
be excited in UHR layer, hence, a second mode conversion takes place from X-mode into an electron Bernstein
wave. Electron Bernstein wave heating and current drive
is an attractive candidate for sustaining a steady state ST
plasma [5, 10, 11]. Accurate understanding of mode conversion processes in fusion devices is important both forthe interpretation of diagnostic data and quantitative analysis of ECH and/or ECCD.
Usual analysis of propagation and absorption of EC
waves employs ray tracing method based on geometrical optics for waves with short wavelength [12, 13]. In a
plasma with high density or low magnetic field, however,the presence of cutoﬀ layer may prevent the waves from
penetrating into the central part from the low field side.
Therefore full wave analysis in which Maxwell’s equation
is solved as a boundary-value problem is required for selfconsistent description of the RF waves. In a hot plasma,
however, we have to include kinetic eﬀects, especially finite Larmour radius (FLR) eﬀects, in the full wave analysis. Several schemes have been developed to describe theFLR eﬀects and implemented in the numerical codes. One
is to expand the FLR eﬀects in the powers of k⊥ ρ wherek⊥ is the wave number perpendicular to the magnetic fieldand ρ is the Larmor radius. This analysis is limited to|k⊥ ρ| < 1 and is usually limited up to the second order, e.g.
the TORIC code [14]. Therefore this scheme is not applicable to the situation where k⊥ ρ  1. The second scheme
is to apply the Fourier analysis to inhomogeneous plasmas
and sum up all the Fourier components of the dielectric ten-2403070-1

c 2016 The Japan Society of Plasma
Science and Nuclear Fusion Research

Plasma and Fusion Research: Regular Articles

Volume 11, 2403070 (2016)

sor. This scheme has been successfully applied to the ion
cyclotron heating using the AORSA code [15]. Since all
Fourier components are coupled with each other, a largesize matrix equation with a dense coeﬃcient matrix has tobe solved. Therefore it is prohibitive to apply this scheme
to EC waves with short wave-length in tokamak plasmas.
We have developed a third scheme to describe the dielectric tensor as an integral operator and solve the integrodiﬀerential Maxwell’s equation [16] by the finite element
method (FEM). Preliminary analysis of this scheme wasapplied to the IC waves in the presence of energetic particle
scheme [17]. Merit of the integro-diﬀerential approach isthat, since the wave number is not used in this approach, it
is applicable to the cases where k⊥ ρ > 1 and that the interaction with wave electric field is spatially localized within
gyro radius of particles, and therefore the coeﬃcient matrix is sparse and less computational resource is requird.
In the present article, the full wave analysis of the electron Bernstein waves in tokamak plasmas is reported forthe first time. Though the analysis is one-dimensional andsmall-size plasma at present, it seems feasible to extend to
two-dimensional and large-size plasma.
In this paper, we introduce one-dimensional kinetic
full wave analysis using the integro-diﬀerential approachin section 2, apply it to the O-X-B mode conversion in
tokamak configuration in section 3, and summary and discussion are given in section 4.

2. Theoretical Formulation of Finite
Larmor Radius Eﬀects
Let us consider the particle trajectory in the presenceof static magnetic field B0 (x) in z-direction at time t, i.e.,v ≡ v(t) = (v⊥ cos φ, v⊥ sin φ), v ), we can write for timet = t − τ as x ≡ x(t ) = x(t) + (v⊥ /Ω)(sin φ − sin(φ + Ωτ))and v ≡ v(t ) = (v⊥ cos(φ + Ωτ), v⊥ sin(φ + Ωτ), v ),where Ω = qB0 /m is the cyclotron frequency, q and m
are the electric charge and mass of the particle, φ is theinitial gyro phase, and v⊥ and v are components of theelectron velocity perpendicular and parallel to the magnetic field, respectively. Considering the electric field perturbation: E(r, t) = E(x) exp ( i k z − i ωt), the linearized
Vlasov equation leads to the induced current expressed as3/2 
1/2T ⊥ (x0 )mn0 (x0 )q2j(x) =T ⊥ (x0 ) 2πT ⊥ (x0 )T  (x0 ) ∞ ∞ 2π  ∞
×dv
dv⊥
dφ
dτ
−∞
000⎧
⎫
⎪
⎪
mv2
⎪
⎪
mv2⊥
⎨
⎬
× exp ⎪−+ i (ω − k v )τ⎪−⎪⎪⎩ 2T ⊥ (x0 ) 2T  (x0 )⎭⎞⎛⎜⎜⎜v⊥ cos φ⎟⎟⎟ ⎟⎜(3)× ⎜⎜⎜⎜⎜ v⊥ sin φ ⎟⎟⎟⎟⎟ G x Gy Gz · E(x ),⎠⎝v
wheredv qv

dτ

where x0 is the guiding center position.
anisotropic Maxwellian distribution function(1)
Using the
3/2 
1/2mT ⊥ (x0 )f0 (x0 , v) = n0 (x0 )2πT ⊥ (x0 )T  (x0 )⎫⎧2 ⎪2⎪mv ⎪⎪mv⊥⎬⎨−× exp ⎪,−⎪⎪⎭⎩ 2T ⊥ (x0 ) 2T  (x0 ) ⎪
T ⊥ k v 2Gx = 1 − 1 −v⊥ cos(φ + Ωτ)T ω
T ⊥ k v 2Gy = 1 − 1 −v⊥ sin(φ + Ωτ)T ω∂T⊥ 2T ⊥ v2⊥ v
Gz =cos(φ + Ωτ)  .
v v − i 1 −T ⊥T ω∂xHere, the transformation from the velocity space variables(v⊥ , φ) to the particle position x and the guiding center position x0 is of prime importance. The relation between them is established through the unperturbed orbit,
x = x0 − (v⊥ /Ω) sin φ and x = x0 − (v⊥ /Ω) sin(φ + Ωτ).
The Jacobian of this transformation is given by J ≡−Ω2 /v⊥ sin Ωτ. The integrand in eq. (3) is a periodic function w.r.t. Ωτ which can be expressed as Fourier series
and easily integrated over τ. The integration over v canbe expressed by the√ plasma dispersion function Z (η) withη ≡ (ω − nΩ) /( 2k vT  ) where vT = (T/m)1/2 . Finally,
the expression for induced current becomes ∞σ (x, x ) · E(x ),j(x) =dx ↔
(4)−∞and the explicit form of conductivity tensor ↔
σ (x, x ) is↔
given in Appendix. The dielectric tensor (x , x) is thesum of the conductivity tensor of species s and the contribution from the displacement current,↔
↔(x , x) =δ(x − x ) + i
σ (x, x ).
ω0 s tq i (ω−k v )τ
em−∞ ∂ f0 (x0 , v )× E(x ) + v × B(x ) ·,∂v

j(x) =

the induced current density is reduced to(5)Using this dielectric tensor, Maxwell’s equation is expressed asω2 ∞  ↔
∇ × ∇ × E(x) − 2dx (x, x ) · E(x )c −∞= i ωμ0 jext (x).
(6)

3. 1D Kinetic Full Wave Analysis ofO-X-B Mode Conversion(2)

We have solved eq. (6) in a slab model numericallyby using the finite element method and analyzed the mode
conversion and power absorption. The plasma parameters2403070-2

Plasma and Fusion Research: Regular Articles

Volume 11, 2403070 (2016)

Fig. 1 (a) Wave number k x as a function of major radius x = R − R0 obtained from the dispersion relation. The red and blue marks indicate
the real and imaginary parts of the wave number, respectively. The marks in the upper half plane corresponds to propagating waves,k2x > 0, and the lower half plane evanescent waves, k2x < 0. (b) Profiles of the wave electric field and the power deposition along
x major radius for T (0) = 500 eV, kz = 32 m−1 . The red line indicates the real part of the wave electric field, the blue line theimaginary part, and the black line the absolute value, respectively.

of the small-size spherical tokamak (LATE) are considered. Major radius of the device R0 = 0.22 m, minor radiusa = 0.16 m, central magnetic field is B0 = 0.08 T, central
electron density is 1017 m−3 and the frequency of the wave
excited by antenna current is 2.45 GHz. The temperature
profile is parabolic and the density profile is the square rootof the parabolic. In the following calculation, 1000 meshpoints are used.
Figure 1 (a) shows the wave number k x as a functionof major radius x = R − R0 obtained from the dispersion relation for the optimum parallel wave number kz = 32 m−1 .
The vertical axis indicates |k x | for propagating modes in
the upper half plane and −|ik x | for evanescent modes in thelower half plane. For this value of kz , both the plasma cutoﬀ (PC: ω2 = ω2p ) of the O-mode and the left-hand-cutoﬀ(LC) of the X-mode are very close to each other aroundx = 0.085 m. The O-mode is mode-converted to the Xmode at PC & LC and the X-mode is mode converted to theEBW near the UHR after changing the direction of propagation (mode conversion from the fast X-mode to the slow
X-mode). The EBW propagates inwards with increasing
wave number towards the high field side, and is absorbed

at the ECR.
The wave electric filed structure and the power deposition profile calculated by the TASK/W1 code for theparameters of Fig. 1 (a) are shown in Fig. 1 (b). For input
power of 1 W, the unit of the wave electric field is V/m, andthat of absorbed power density is W/m3 . The wave is excited as an O-mode by the antenna current in the z direction
at R−R0 = 0.17 m. The parallel component Ez indicates thesum of the O-mode component and the electrostatic component − i kz Φ of the wave. Most of the O-mode is modeconverted to the X-mode component Ey and the other part
penetrate over the evanescent layer to the high field side.
Though a part of the X-mode returns to the plasma surface
x = 0.16 m, most of the wave energy is converted to theelectrostatic Bernstein mode near the UHR, x = 0.12 m.
The electrostatic component is dominant in the x component of the wave electric field E x , since the wave lengthof the Bernstein is very short and the wave number k x is
much higher than those of the X- and O-mode. The Bernstein wave propagates towards the high density region andabsorbed near the fundamental ECR around of x = −0.1 m.
The absorbed power density pabs (x) is simply calculated by

2403070-3

Plasma and Fusion Research: Regular Articles

Volume 11, 2403070 (2016)

Pabs = j∗ · E and oscillates owing to the interference between the EBW and the X- mode. This comes from theambiguity between the kinetic power absorption and thegradient of kinetic power flux. More rigorous definition ofPabs is necessary for positive power absorption density.
The injection angle dependence of the mode conversion eﬃciency is examined by the kz dependence of the antenna loading resistance. Preliminary results suggest thatthe antenna loading resistance has a maximum at the optimum angle, kz = 32 m−1 , and decreases to about 30% forkz = 24 m−1 and kz = 40 m−1 .

4. Summary
We have made kinetic full wave analysis of electron
cyclotron waves for O-X-B mode conversion by 1D hot
plasma model by employing the TASK/W1 code. Suchanalysis using the integral form of dielectric tensor derived
by following the unperturbed particle orbit successfully describes the O-X-B mode conversion near the UHR layer
and absorption near the ECR layer. More systematic analyses of parameter dependence and quantitative comparison
with analytical estimates will be reported in a separate paper.

Acknowledgment
This work is supported by Grant-in-Aid for Research
Fellows with JSPS KAKENHI Grant Number 26-04380
from JSPS, Japan.

Hyz = A2 (X − Y)Fn(1)√2vT  ηHzz =A2 Fn(1) ,vT ⊥withX≡Ω
x + x
x0 −,vT ⊥2Y≡Ω(x − x ).
2vT ⊥The kernel function Fn(i) is defined by π1Fn(i) (X, Y) ≡ 2dθ
2π 0
Y2X2−× exp −f (i) (θ),
1 + cos θ 1 − cos θ nwith⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨(i)fn (θ) = ⎪⎪⎪⎪⎪⎪
⎪
⎪
⎪⎪⎪⎪⎪
⎩

cos nθ
sin θsin nθ
sin nθ
sin2 θcos θ sin nθ
sin2 θ

(i = 1)(i = 2)(i = 3)(i = 4).

Finally the coeﬃcients A1 and A2 are given byωT ⊥ Z  (η),A1 ≡ √Z(η) + 1 −T22k vT T⊥ω
ΩT⊥
A2 ≡+nZ  (η)·
1−2k vT ⊥ T ωT

Appendix
The integral form of the conductivity tensor is writtenas↔
σ (x, x ) =2 
∞i n0 q2 Ω↔
H (x, x , x0 ),
dx0mω
v
T⊥
−∞n=−∞

 ∞(A.1)↔
where the components of the tensor H are given by
H xx = −nA1 Fn(2)Hyx = iA1 (X − Y) (X − Y)Fn(3) − (X + Y)Fn(4)
Hzx = −iA2 (X − Y)Fn(3) − (X + Y)Fn(4)H xy = −iA1 (X + Y) (X + Y)Fn(3) − (X − Y)Fn(4)Hyy = −A1 (X + Y)(X − Y)Fn(1)Hzy = A2 (X + Y)Fn(1)H xz = iA2 (X + Y)Fn(3) − (X − Y)Fn(4)[1] N.J. Fisch, Rev. Mod. Phys. 59, 175 (1987).
[2] R. Prater, J. Fusion Energy 9, 19 (1990).
[3] R. Prater, Phys. Plasmas 11, 2349 (2004).
[4] B. Lloyd, Plasma Phys. Control. Fusion 40, A119 (1998).
[5] J. Urban et al., Nucl. Fusion 51, 083050 (2011).
[6] J. Preinhaelter and V. Kopecky, J. Plasma Phys. 10, 1 (1973).
[7] H.P. Laqua et al., Phys. Rev. Lett. 78, 3467 (1997).
[8] H.P. Laqua, Plasma Phys. Control. Fusion 49, R1 (2007).
[9] S. Shiraiwa et al., Phys. Rev. Lett. 96, 185003 (2006).
[10] H. Idei et al., J. Plasma Fusion Res. SERIES 8, 1104 (2009).
[11] H. Igami et al., Nucl. Fusion 49, 115005 (2009).
[12] T. Maekawa et al., J. Phy. Soc. Jpn. 48, 247 (1980).
[13] A. Fukuyama, Fusion Eng. Des. 53, 71 (2001).
[14] M. Brambilla, Plasma Phys, Control. Fusion 41, 1 (1999).
[15] E.F. Jaeger et al., Phys. Plasmas 9, 1873 (2002).
[16] O. Sauter and J. Vaclavik, Nucl. Fusion 32, 8 (1992).
[17] A. Fukuyama et al., Proc. 21st Fusion Energy Conf.,
Chengdu, China, 2006, TH0P6-10, IAEA (2006).

2403070-4