Plasma and Fusion Research: Regular Articles

Volume 19, 1403026 (2024)

Two-Dimensional Axisymmetric Finite Element Simulation ofLower Hybrid Wave with an Iterative Scheme
Fumiya ADACHI, Naoto TSUJII, Akira EJIRI, Kouji SHINOHARA, Seowon JANG, Yi PENG,
Kotaro IWASAKI, Yu-Ting LIN, Zhengnan JIANG, Yiming TIAN, Yangguang JIANG,
Shengyu WANG, Yijin XIONG and Masaru YOSHIDA
University of Tokyo, Kashiwa 277-0882, Japan(Received 27 October 2023 / Accepted 25 June 2024)
Simulating waves in hot plasmas in configuration space is a diﬃcult problem because of the non-local property of the plasma response, which makes the wave equation an integro-diﬀerential one. In this research, we
conducted an axisymmetric hot plasma full wave simulation at the lower hybrid frequency range with the finite
element method. Kinetic eﬀects were introduced in the direction parallel to the magnetic field. We implemented
the code so that it can handle an arbitrary velocity distribution function to introduce electron kinetics. An iterative method was utilized to introduce the non-local hot plasma contribution, which is more memory eﬃcient thandirect solving. The hot plasma perturbed current density was iteratively calculated in our scheme. For the present
simulation, we used the plasma equilibrium obtained by the experiment with the TST-2 spherical tokamak, which
is located at the University of Tokyo. The simulation with the realistic profiles successfully converged. The
electron heating power deposition profile was estimated for the obtained electric field solutions.
c 2024 The Japan Society of Plasma Science and Nuclear Fusion Research⃝

Keywords: lower hybrid wave, full wave simulation, finite element method
DOI: 10.1585/pfr.19.1403026

1. Introduction
Accurate analyses of waves in plasmas need a wave
theory that takes into account the kinetic response of theplasma, namely, the hot plasma model. The diﬃcult pointof the hot plasma full wave simulations is that the plasma
response becomes non-local, expressed as an integral operator. This makes the wave equation to be solved an integrodiﬀerential equation in configuration space. Because ofthis characteristic, many of the previous full wave calculations were based on the spectral representation [1, 2].
The plasma response becomes an algebraic product inwavenumber space, making the implementation straightforward. This natural choice of simulation domain, however, requires the inversion of a large and dense matrix, andthe problem becomes computationally expensive. Also,
owing to the use of a global Fourier basis, it is ineﬃcient
to model complex geometries such as the scrape oﬀ layer
and the antenna structure. In real experiments, interactionsof waves with plasmas in the scrape oﬀ layer and/or withlaunchers could lead to deteriorated eﬃciency of heating
and current drive. The interactions need to be simulated
accurately. Full wave simulations in configuration spacewith the finite element method (FEM) can eﬃciently handle complicated structures with flexible meshing, such as
using finer meshing where more meshing is needed. The
FEM enables integrated modelling of the core plasma andedge regions.

Lower hybrid (LH) current drive has been conductedin many tokamak experiments [3, 4]. Traditionally, ray
tracing has been used to analyze waves in the LH frequency range. The technique solves the geometrical optics
equations assuming the Wenzel-Kramers-Brillouin (WKB)
approximation. It succeeded, for example, in explaining
experimentally obtained current drive eﬃciency [5]. However, coupling of LH waves from the antenna to the plasma
cannot be treated within the framework because the LH
wave has a low-density cutoﬀ at ω = ωpe , where ω/2π isthe frequency of the LH wave and ωpe /2π is the electron
plasma frequency. Ray tracing must start somewhere inside the plasma where the LH waves are not evanescent. In
addition, when the LH waves approach the cutoﬀ layer, thewavelength and the plasma scale length become comparable and the WKB approximation may break down. This
could be a problem in the multi-pass damping regime. Fullwave simulations can overcome these diﬃculties.
Previous FEM full wave simulations solved theintegro-diﬀerential wave equation with an iterative scheme[6–8]. An iterative scheme consumes more time than adirect solver, but can work with less computational resources. In Ref. [6], two-dimensional finite element fullwave simulations of the LH waves were implemented. A
possible problem is that their scheme included division by
the electric field to calculate an eﬀective dielectric. It may
be inaccurate at interferences where the electric field becomes almost zero. In this work, we present a diﬀerent

author’s e-mail: adachi@fusion.k.u-tokyo.ac.jp

1403026-1

c 2024 The Japan Society of Plasma⃝
Science and Nuclear Fusion Research

Plasma and Fusion Research: Regular Articles

Volume 19, 1403026 (2024)

iterative scheme to simulate a two-dimensional LH wave
model. In our scheme, the hot plasma perturbed current
density that is consistent with the global wave simulation
is iteratively calculated.
The organization of this paper is as follows. Section 2
briefly describes the wave physics and the implementationof our code. Section 3 shows the results from the homogeneous hot plasma simulation of LH waves and verifies thepresent scheme. In Sec. 4, we simulate LH waves in a hot
plasma with a realistic plasma equilibrium and show theresults. Finally, in Sec. 5, we summarize our work.

2. Implementation of LH Wave Model
Assuming a time dependence ∝ eiωt , high frequency
waves are described by∇ × (∇ × E) −ω2
E = −iωµ0 j,c2(1)

where E is the electric field, j is the current density, c is the
speed of light in vacuum and µ0 is the vacuum permeability. In wave analysis, we usually assume a linear plasma
response to waves:j(r) = σ̂ · E.

(2)

Here, σ̂ is the plasma conductivity. With Eq. (2), Eq. (1)becomes the Helmholtz equation:∇ × (∇ × E) −

ω2
K̂ · E = 0,
c2

(3)

where
K̂ = 1 −

iσ̂
,
ϵ0 ω

(4)

K̂ is the plasma dielectric and ϵ0 is the vacuum permittivity. The plasma dielectric or the plasma conductivity determines wave propagation and damping. In general, theplasma dielectric and the plasma conductivity are integral
operators. Therefore, the wave equations Eqs. (1) and (3)become integro-diﬀerential equations. However, when we
consider a uniform plasma or assume local homogeneity,
the plasma dielectric and the plasma conductivity become
algebraic multiplication in the spectral domain.
The major linear damping mechanism for the LH
waves is electron Landau damping. The corresponding
term of the conductivity can be obtained as below by
following the derivation of [9]. Here, we consider onedimensional Vlasov equation of species s:∂ fs∂ fs qs E ∂ fs
+ vz
+
= 0,
∂t
∂z
m s ∂vz

Equation (6) can be solved by integrating along the unperturbed particle orbit:∫ t∂ f s0 ′qs(z , vz ), (7)dt′ E(z′ , t′ )f s1 (z, vz , t) = −m s −∞∂vz
where z − z′ = vz (t − t′ ). By taking the velocity moment,
the perturbed current density becomes
(
)
∫ ∞
∫
q s 2 ∞ z − z′ ∂ f s0 iωτ
j s1 (z) =
dz′ −
dτ 2
e
E(z′ ).
m
∂v
τ
s
z
−∞
0
(8)
To obtain a simplified expression of the conductivity operator in configuration space, we have changed the variablesof integration: (t′ , vz ) → (t′ , z′ ) and set τ = t − t′ .
When we assume local homogeneity f s0 (z′ , vz ) ≈f s0 (z, vz ), the conductivity kernel is∫
q s 2 ∞ z − z′ ∂ f s0 iωτ
′
dτ 2
e .
(9)
σhot,s (z − z ) = −
ms 0
τ ∂vz
The real part of Eq. (9) corresponds to electron Landau
damping, and the imaginary part introduces a kinetic correction in wave propagation. Note that in Eq. (8), the current density j s1 is obtained as a convolution of the conductivity kernel and the electric field. With ω/2π = 200 MHz,an electron density of 3 × 1017 m−3 , and a velocity distribution function shown in Fig. 1, the conductivity kernel becomes as shown in Fig. 2. The velocity distribution function has two characteristic temperatures: the bulk temperature T b = 15 eV and the tail temperature T t = 20 keV. For
vz < 1.5vth,b , the distribution follows a Maxwellian distribution with temperature T b , and for vz > 1.5vth,b , the distribution follows a Maxwellian distribution with temperature
T t . Here, vth,b is the thermal velocity of bulk electrons, de√fined as vth,b = 2T b /me .

(5)where f s (z, vz , t) is a velocity distribution function, q s is acharge, m s is a mass. Considering a perturbation f1s aroundthe equilibrium distribution f0s , Eq. (5) can be linearized to()∂∂
q s E ∂ f s0
d f s1=+ vz
f s1 = −
.
(6)
dt
∂t
∂z
m s ∂vz

Fig. 1 Example of the electron velocity distribution function
with the electron density of 3 × 1017 m−3 . It switchesfrom the Maxwellian with the bulk temperature of 15 eV
to the Maxwellian with the tail temperature of 20 keV.
The switching velocity is 1.5 times the bulk thermal velocity in this case.

1403026-2

Plasma and Fusion Research: Regular Articles

Volume 19, 1403026 (2024)
where k∥0 is a prescribed parallel wavenumber andσhot,e (k∥ ) is the Fourier transform of σhot,e (z). jNL is the
fully non-local expression of the current density and jA is
the approximate expression of the current density with theprescribed wave number. Also, we define δj as the diﬀerence between those two quantities:δj = jNL − jA .

Fig. 2 Example of the conductivity kernel expressed in Eq. (9).
The wave frequency of ω/2π = 200 MHz, and the electron velocity distribution function as shown in Fig. 1 areassumed.

Our simulation used the plasma dielectric which introduced electron Landau damping:K̂ = Kcold,⊥ + K̂hot,∥


 S
iD 0= −iD S 0 + ẑẑ(P + χ̂ELD ),


0
0 0

(10)

where S , D and P are the Stix parameters as defined in [10],andiσ̂ELDχ̂ELD = −
,
(11)
ϵ0 ω
σ̂ELD = ℜσ̂hot,e .
(12)Equation (10) is written in the magnetic coordinate system
where the magnetic field is in the ẑ direction. In the cylindrical coordinates (R, ϕ, Z), the plasma dielectric Kc can beobtained by rotating the coordinates as
Kc = U−1 · K · U,
where

√

2
 1 − bR
U = 
0


 b
R

(13)
−bR bϕ√1 − bR 2bZ√
1 − bR 2bϕ

−bR bZ √
1 − bR 2 −bϕ  ,
√
1 − bR 2 

bZ

(14)

bR , bϕ , bZ are the three components of the normalized magnetic field in the cylindrical coordinates:
B= bR R̂ + bϕ ϕ̂ + bZ Ẑ,
|B|

(15)

and B is the magnetic field.
Here we introduce two expressions for the perturbed
current density:∫
jNL = ẑ dz′ σELD (z − z′ )Ez (z′ ),(16)jA = ẑσELD (k∥0 )Ez ,

(17)

(18)

Since Eqs. (1) and (3) are integro-diﬀerential equations, wesolved the wave equation iteratively as follows:
 (
) ω2




∇× ∇×E(N) − 2 (Kcold,⊥ + Khot,∥ (k∥0 ))·E(N)



c

(N−1)


=
−iωµ
I
0


(
)



I(N−1) = A δj(N−1) , δj(N−2) , . . . , δj(0) ,(19)
where N is the iteration index, Khot,∥ (k∥ ) is the Fouriertransform of the kernel Khot,∥ (z), and A symbolizes theoperation of Anderson acceleration [11]. In order to improve convergence, the non-local plasma eﬀect was split
into the current density and the plasma dielectric. Only
the diﬀerence between the non-local expression and theapproximate expression was iterated. In addition, we utilized Anderson acceleration, one of the vector extrapolation methods. In our simulation, we used all the previously obtained solutions for Anderson acceleration to prepare the next input I at every iteration. The RF module of
COMSOL Multiphysics [12] was used to solve the wave
equation in Eq. (19) with the FEM. By driving the RF
module and evaluating the perturbed current density with
MATLAB [13], we iteratively calculated wave propagationand damping.

3. Homogeneous Plasma Simulations
and Results
First, we simulated LH waves in a homogeneous
plasma described in Fig. 3 to verify the iterative scheme
explained in Sec. 2. The electron density was 1 ×1017 m−3 , and the electron velocity distribution function
was Maxwellian with a temperature of 500 eV. The magnetic field with a magnitude of 0.04 T was tilted for 92◦
from the x̂ direction. This is to simulate the tilted magnetic field in a tokamak because of the poloidal field. The
length of the simulation domain in the ŷ direction was thesame as the wavelength of the LH waves in the ŷ direction, and the parallel wavenumber was scanned by changing the wavenumber in the ŷ direction ky at the boundary.
We added an artificial damping region to avoid the reflection of the wave at the edge.
The waves were excited with the electric field boundary condition; the ŷ component of the electric field was
imposed as Ey = Eant sin(ky y), where Eant = 1 V/m. The
wave frequency was 200 MHz. The imposed wavenumber in the ŷ direction ky was used as a prescribed parallel wavenumber k∥0 . In this tilted magnetic field geome-

1403026-3

Plasma and Fusion Research: Regular Articles

Volume 19, 1403026 (2024)
the wave damping became stronger for the converged solution than the initial solution. Also, as ky increased in
the chosen parameter range, the wave damping became
stronger. This agreed with the characteristic of Landau
damping. The converged solutions were almost the sameas the dispersion relation with wavenumber upshift. The
dispersion relation with the upshift was obtained by solving the dispersion relation withBy
Bx
+ ky .
(20)
B
B
The relation without the upshift was obtained in the sameway, but by neglecting the first term on the right-hand sideof Eq. (20).
We can roughly evaluate wave damping in a hot
plasma by approximating a wavenumber in the direction
parallel to the magnetic field with a toroidal wavenumber(ky in this model) [14]. The parallel wavenumber, however, deviates from the toroidal wavenumber because of
the poloidal magnetic field. It is necessary to introduce thisdiﬀerence into a model for quantitatively accurate simulations. With the present iterative scheme, wave damping
can be calculated with more accuracy by using the integral
form of the hot plasma perturbed current density.
k∥ = k x

Fig. 3 Settings for the homogeneous plasma simulation.

4. Axisymmetric Full Wave Simulation and Results
Fig. 4 Residual error for the homogeneous plasma simulation.

Fig. 5 Imaginary part of the wavenumber in the x̂ direction. The
values for the initial and converged solutions and the analytical values are shown.

try, the parallel wavenumber of the analytical solution islarger than ky , and the wavenumber correction was iteratively added to the solution.
Figure 4 shows the residual error for the homogeneous
plasma simulation. The residual error for all the simulations became suﬃciently small successfully. Figure 5
shows the imaginary part of the wavenumber in the x̂ direction. As expected, since the wavenumber upshift occurred,

TST-2 is a spherical tokamak located at the University of Tokyo. Its major radius is 0.36 m, its minor radius
is 0.23 m, and its toroidal field strength is < 0.3 T. Threecapacitively coupled combline antennas are installed at theouter-midplane, at the top, and at the outer-oﬀ-midplane[15–17]. The parallel refractive index is 5.5 for the outermidplane launch antenna, 4.9 for the top launch antenna,and 13 for the outer-oﬀ-midplane launch antenna, respectively. A source of 100 kW is available for each antenna
operating at a frequency of 200 MHz. Typically, the discharge duration is about 80 ms. The electron density andthe electron temperature profiles are measured with Thomson scattering diagnostics, and the fitted curves at the flat
top of a discharge were used in the present simulation.
We conducted a two-dimensional axisymmetric fullwave simulation on a poloidal plane. Toroidal mode number, which is a conserved quantity in an axisymmetric system, was nϕ = 20 and the wave frequency was ω/2π =200 MHz. The bulk electron temperature and the electron density were assumed to be flux functions. The magnetic flux surfaces were obtained by EFIT [18] reconstruction and is shown in Fig. 6. The bulk electron temperature and the electron density profiles used in the simulation are shown in Fig. 7. These are the equilibrium profiles obtained in the experiment with the TST-2 sphericaltokamak. We used the solutions from a Fokker-Plancksolver CQL3D [19] coupled with ray tracing code GENRAY [20] for the electron velocity distribution function.
Figure 8 shows the used electron velocity distribution func-

1403026-4

Plasma and Fusion Research: Regular Articles

Volume 19, 1403026 (2024)

Fig. 8 Electron velocity distribution function used in the axisymmetric simulation.
This is obtained from the
GENRAY-CQL3D simulation result.

Fig. 6 Normalized toroidal flux obtained by EFIT. The simulated area is encircled by the black solid lines.

Fig. 9 Residual error for the axisymmetric simulation.

Fig. 7 Bulk electron density and electron temperature profiles of
the LH start-up plasma in TST-2.

tion. This is obtained by integrating the GENRAY-CQL3D
simulation result over the perpendicular velocity. Outside
the last closed flux surface (the white region in Fig. 6), thebulk electron temperature and the electron density were assumed to be constant and the same as that on the last closed
flux surface. The waves were excited with the electric fieldboundary condition; the toroidal component of the electricfield was imposed as follows on the outer edge around the
midplane.
Eϕ = 1 V/m (−0.5L < z < 0.5L),(21)

with L = 0.15 m. The other walls were assumed to be
perfect electric conductors. We initialized the iteration
with I(−1) = 0, that is, by considering only the approximate analytic damping. To calculate the approximate current density, Eq. (17), we prescribed toroidal wavenumber
ktor = nϕ /R as an approximation for the parallel wavenumber. The current density component corresponding to thediﬀerence between σhot,e (k∥ ) and σhot,e (ktor ) was iteratively
added to the model.
Figure 9 shows the residual error for the axisymmetricsimulation. The relative change in the total power deposition is shown in Fig. 10. Here P(N) is the electron heating
power density ℜ(jNL · E∗ )/2 integrated over the simulation
domain at the Nth iteration. Although the error of the total
power deposition was fluctuating, it still stayed at the orderof 10−4 after about 15 iterations. The converged solution
is shown in Fig. 11. Figure 12 shows the profile of theheating power density averaged over the toroidal flux surfaces at the initial iteration and the 20th iteration, and theresult from a ray tracing code which only considers electron Landau damping. In the simulated plasma, the parallel

1403026-5

Plasma and Fusion Research: Regular Articles

Volume 19, 1403026 (2024)

Fig. 10 Relative change of the deposited power. The change
was less than 10−3 after 17 iterations.

wavenumber of the waves experienced a downshift due tothe poloidal magnetic field [21]. The actual wave absorption must be weaker than the approximate analytic absorption (= initial absorption). Therefore, as the iteration wenton, the region where the wave was damped moved outward. Note that the waves were launched from the low fieldside and damped near the inner walls with a single pass inthe present simulation. The converged solution from thefull wave simulation showed absorption near the plasmaedge, and the characteristic matched that of the ray tracing
result. However, they showed a slightly diﬀerent locationof absorption. Full wave eﬀects such as diﬀraction are notincluded in ray tracing. Also, it is non-trivial how to set
the initial wavenumber for each ray in ray tracing, whichaﬀects the location of absorption. The discrepancy may
come from these reasons.

5. SummaryFig. 11 Converged electric field solution. The parallel component is shown.

Fig. 12 Flux surface averaged electron heating power density
profile.

We have developed a new LH full wave simulation
code with the FEM scheme, utilizing the commercial FEM
solver. The new code is able to handle an arbitrary electron
velocity distribution function. The non-local hot plasma
contribution was implemented with an iterative scheme.
The electron Landau damping was introduced as a nonlocal hot plasma response. In order to improve convergence, we split the hot plasma eﬀect into the plasma dielectric and the current density, and made use of Anderson acceleration. With the newly developed code, firstly
we have conducted a homogeneous plasma simulation andverified the code. Secondly, we have conducted an axisymmetric simulation of the LH waves with the realistic electron bulk density and temperature profiles. The simulation
successfully converged. The electron heating power profile was estimated for the converged electric field and theresult agreed well with that of the ray tracing code.
Interactions of waves with plasmas in the scrape oﬀ
layer and/or with launchers could be problematic for heating and current drive. It is necessary to consider the interactions in real experiments and future reactors. The previous codes using spectral representation cannot eﬃciently
describe physics in such complex geometries. The FEM
enables integrated modelling of the antenna structure andthe hot core plasma, but the integral formulation of plasma
response must be employed. We introduced the integral
form perturbed current density in the present FEM model.
In the future extensions, the integrated modelling with our
FEM model will be tested. Also, the present model neglects the kinetic propagation correction and handles only
the real part of the conductivity kernel. The full treatmentof the kinetic expressions remains as a future work.
[1] M. Brambilla and T. Krücken, Nucl. Fusion 28, 1813
(1988).
[2] E.F. Jaeger, L.A. Berry, E. D’Azvedo, D.B. Batchelor et al.,Phys. Plasmas 8, 1573 (2001).

1403026-6

Plasma and Fusion Research: Regular Articles

Volume 19, 1403026 (2024)[3] S. Shiraiwa, S. Ide, S. Itoh, O. Mitarai et al., Phys. Rev.
Lett. 92, 035001 (2004).
[4] V.V. Dyachenko, O.N. Shcherbinin, E.Z. Gusakov, V.K.
Gusev et al., Nucl. Fusion 55, 113001 (2015).
[5] P.T. Bonoli, J. Ko, R. Parker, A.E. Schmidt et al., Phys.
Plasmas 15, 0561117 (2008).
[6] O. Meneghini, S. Shiraiwa and R. Parker, Phys. Plasmas
16, 090701 (2009).
[7] D. Green and L. Berry, Comput. Phys. Commun. 185, 736(2014).
[8] P. Vallejos, T. Hellsten and T. Johnson, J. Phys. Conf. Ser.
1125, 012020 (2018).
[9] J.E. Drummond, R.A. Gerwin and B.G. Springer, J. Nucl.
Energy, Part C Plasma Phys. 2, 98 (1961).
[10] T. Stix, Waves in Plasmas (American Institute of Physics,
New York, 1992).
[11] H. Walker and P. Ni, SIAM J. Numer. Anal. 49, 1715(2011).
[12] COMSOL Multiphysics, www.comsol.com, COMSOLAB, Stockholm, Sweden.
[13] MATLAB, www.mathworks.com, The MathWorks Inc.,
Natick, MA, USA.
[14] C. Lau, L.A. Berry, E.F. Jaeger and N. Bertelli, Plasma
Phys. Control. Fusion 61, 045008 (2019).
[15] T. Shinya, Y. Takase, S. Yajima, C. Moeller et al., Nucl.
Fusion 57, 036006 (2017).
[16] S. Yajima, Y. Takase, Y. Tajiri, Y. Takei et al., Nucl. Fusion
59, 066004 (2019).
[17] Y. Ko, N. Tsujii, A. Ejiri, O. Watanabe et al., Nucl. Fusion
63, 126015 (2023).
[18] L.L. Lao, H. St. John, R.D. Stambaugh, A.G. Kellman etal., Nucl. Fusion 25, 1611 (1985).
[19] R.W. Harvey and M.G. McCoy, Proc. IAEA TCM on Advances in Sim. and Modelling of Thermonuclear Plasmas,
489 (1992).
[20] A.P. Smirnov, R.W. Harvey and K. Kupfer, Bull. Am. Phys.
Soc. 39, 1626 (1994).
[21] P.T. Bonoli and E. Ott, Phys. Fluids 25, 359 (1982).

1403026-7