第 17 回
核燃焼プラズマ
統合コード研究会
17th Burning Plasma Simulation
Initiative (BPSI) Meeting
日時：2019 年 12 月 5 日(木)－6 日(金)場所：九州大学筑紫キャンパス
応用力学研究所 2 階大会議室(Ver. 1)

合同会合第 17 回核燃焼プラズマ統合コード研究会
17th Burning Plasma Simulation Initiative (BPSI) Meetingおよび
核融合エネルギーフォーラムサブクラスター会合
シミュレーションクラスター 炉心プラズマモデリングサブクラスター
プラズマ物理クラスター 定常運転・制御サブクラスター(Ver.1)
日時：2019 年 12 月 5 日(木)－6 日(金)場所：九州大学筑紫キャンパス 応用力学研究所 2 階大会議室

12 月 5 日(木)9:15 – 9:25 はじめに
村上 (京大)(座長：村上)9:25 – 9:50 講演 1-1
本多 (量研)Current status of the integrated model GOTRESS+
9:50 – 10:15 講演 1-2
仲田 (核融合研)
Global turbulent transport simulation using a coupled code TRESS+GKV
10:15 – 10:40 講演 1-3
登田 (核融合研)
Reduced model of turbulent transport for helical plasmas including effects ofzonal flows and trapped electrons
10:40 – 10:55 休憩(座長：糟谷)10:55 – 11:20 講演 1-4
小菅 (九大)Interplay among 3d flows in turbulent plasmas
11:20 – 11:35 講演 1-5
青木 (九大)Neutral effects on the structure of minimum enstrophy flows in plasmas
11:35 – 11:55 講演 1-6
山本 (京大)Effects of electron cyclotron heating on the toroidal flow in HSX plasmas
11:55 – 12:20 講演 1-7
村上 (京大)Modeling of supra-thermal electron flux and toroidal torque by ECH in
non-axisymmetric toroidal plasmas
12:20 – 12:30 議論
12:30 – 13:30 昼休み
1

(座長：林)13:30 – 13:55 講演 2-1
相羽 (量研)Impact of plasma rotation and ion diamagnetic drift effects on MHD stability ina QH-mode plasma
13:55 – 14:10 講演 2-2
富松 (九大)Nonlinear dynamics of resistive ballooning instability in PLATO tokamak
plasma14:10 – 14:35 講演 2-3
矢木 (量研)Particle transport model in plasma peripheral region
14:35 – 15:00 講演 2-4
矢本 (量研)Study of multiple impurity seeding effect using SONIC integrated divertor code
for JT-60SA plasma prediction
15:00 – 15:25 講演 2-5
大澤 (九大)Stability of divacancy and interaction with inevitable impurities in metals
15:25 – 15:40 講演 2-6
下村 (九大)Tungsten transport analysis using integrated code TASK
15:40 – 15:55 休憩(座長：本多)15:55 – 16:20 講演 2-7
糟谷 (九大)Analysis of laboratory torus plasma using the integrated simulation scheme
16:20 – 16:45 講演 2-8
沼波 (核融合研)
Development of transport prediction scheme based on first-principle simulationand machine learning
16:45 – 17:10 講演 2-9
成田 (量研)Neural-network-based multi-channel turbulent transport modeling
17:10 – 17:25 講演 2-10
森下 (京大)Data assimilation system based on integrated transport simulation of LHD
plasma17:25 – 17:45 議論
17:45

散会

18:30 – 21:00 懇親会 (はるら 西鉄白木原駅近)2

12 月 6 日(金)9:15 – 9:20 事務連絡
(座長：小菅)9:20 – 9:45 講演 3-1
Chio Zong Cheng (東大)TAEs and Alfven-slow eigenmodes in tokamaks
9:45 – 10:10 講演 3-2
Fukuyama 福山 (京大)Progress in wave heating and current drive using the TASK code
10:10 – 10:35 講演 3-3
Sasaki 佐々木 (九大)Extraction of plasma turbulence dynamics by using dynamical mode
decomposition
10:35 – 10:50 休憩(座長：沼波)10:50 – 11:15 講演 3-4
横山 (核融合研)
Modelling of thermal transport in fusion plasmas based on statistical approach
11:15 – 11:30 講演 3-5
太田 (京大)Evaluation of the momentum dependence of radial diffusion coefficient on
density and heat transport in tokamaks
11:30 – 11:50 議論
11:50 – 13:00 昼休み(座長：藤田)13:00 – 13:05 サブクラスター関係連絡・確認事項
林 (量研)
13:05 – 13:25 講演 4-1
宮戸 (量研)IFERC 計算機シミュレーションセンターの現状報告
13:25 – 13:55 講演 4-2
本多 (量研)ITER に於ける統合コード（IMEG）活動報告
13:55 – 14:35 講演 4-3
若月、鈴木、林 (量研)
ITPA 統合運転シナリオグループ活動報告
14:35 – 14:45 休憩(座長：打田)14:45 – 15:25 講演 4-4
鈴木、林 (量研)JT-60SA の研究活動報告
15:25 – 16:00 サブクラスター活動に関する総合討論と今後の予定
16:00

散会
3

林 (量研)

Current status of the integrated model GOTRESS+
M. Honda, N. Aiba, *H. Seto, E. Narita and N. Hayashi
National Institutes for Quantum and Radiological Science and Technology, Naka,
Ibaraki 311-0193 Japan*National Institutes for Quantum and Radiological Science and Technology, Rokkasho,
Aomori 039-3212 Japan

Recent progress in GOTRESS
Due to the sensitivity of turbulent fluxes to profile gradients, robust predictions arestill arduous with a stiff turbulent transport model incorporated in a conventional transport
code. To overcome this difficulty, we have developed a steady-state transport equation
solver, GOTRESS, which benefits from global optimization techniques such as a genetic
algorithm [1] and the Nelder-Mead method [2] such that they directly find the solution
where a transport flux matches an integrated source [3]. In addition to quick calculationof steady-state profiles, GOTRESS has advantage to make it possible to deal with TGLF
transport model [4], which is one of stiff transport models, in a robust manner on
Multiple-Program Multiple-Data (MPMD) structure framework. Also, GOTRESS has anaffinity for Machine-Learning techniques because it takes advantages of global
optimization techniques that in typical produce a large amount of data, adequate fortraining of an artificial neural network model [5]. This surrogate model mimicking the
behavior of an original model realizes rapid computation.

Novel integrated model GOTRESS+
A novel integrated model GOTRESS+ has been developed, mainly consisting of theequilibrium and current profile solver ACCOME and the neutral beam (NB) heating code
OFMC, other than GOTRESS, for predictions of temperature profiles consistent with anequilibrium and heating profiles. A procedure of GOTRESS+ execution is described in
the following. Figure 1 demonstrates a workflow of GOTRESS+ in the case of high-beta
steady-state scenario in JT-60SA, what is called #5-1: 𝐵𝑇 = 1.72T , 𝐼𝑝 = 2.3MA ,𝑃𝑁𝐵 = 18.85MW and 𝑃𝐸𝐶 = 7MW [6]. The settings of auxiliary heating systems suchas NB heating and electron cyclotron heating (ECH) and the coils that support anequilibrium are described in an input file of ACCOME.

The prescribed density andtemperature profiles are also given as initial guess. After iterative calculations in
ACCOME, a consistent solution is obtained between an equilibrium and a current profile
according to the given kinetic profiles. At the same time, ACCOME also estimates ECH
and EC current drive (ECCD), if needed. The profiles required for OFMC as input aresent to OFMC, and then the profiles of NB heating and NB current drive (NBCD) are
estimated. Now the current density profiles, the equilibrium and heating profiles havebeen computed, which are all requisites for GOTRESS. GOTRESS in turn predicts thetemperature profiles based on them. At this point, the first iteration of GOTRESS+ iscompleted, shown as ``iteration 1'' in fig. 1. The predicted temperature profiles are sent
to ACCOME and replaces the prescribed ones used in the first iteration. The second
iteration then commences. This iterative calculation continues until the profiles are wellconverged. For this case, all the profiles are converged after 6 iterations. The developed
plasma meets the target values: the normalized beta of 4.42 > 4.3 , the energy

Fig. 1: Workflow of GOTRESS+. The predicted profiles are illustrated and each color of a line
corresponds to the number of times of iteration.

confinement improvement of 1.65 > 1.3 and the loop voltage of 0.003~0.
We now compare steady-state profiles of the temperature and the heating power
predicted by GOTRESS+ and TOPICS with the CDBM turbulent transport model
exploited [7]. The boundary conditions were imposed at 𝜌 = 0.8 in both simulations,where 𝜌 is the normalized radial coordinate, so as not to include the pedestal because
the physics governing the formation of the pedestal is different from that in the core region.
Regardless of completely different numerical schemes used, both codes clearly predicted
almost the same temperature profiles. The GOTRESS+ simulation took 6,812 seconds,
around 10 hours faster than the TOPICS simulation. Thanks to the easy-to-use feature,
GOTRESS readily makes it possible to lower the loop voltage much closer to nil bymanipulating heating power carefully.
For more self-consistent predictions, estimation of the pedestal height and width isnecessary. For this purpose, The EPED1 model has widely been used to predict pedestal
height and width [8]. GOTRESS+ is now being extended to include our own EPED1model that uses MARG2D, which is the ideal MHD stability code, replacing the ELITE
code originally used in [8]. In EPED1, temperature profiles are scaled up/down if theinitial ones render a plasma MHD stable/unstable and then an equilibrium is reevaluatedusing ACCOME. This process is repeated until the marginal stability condition is met.
Coupling of GOTRESS+ and EPED1 is in progress.

References[1] P. Charbonneau, Astrophys. J. Suppl. Ser. 101 (1995) 309.
[2] J. A. Nelder and R. Mead, Computer J. 7 (1965) 308.
[3] M. Honda, Comput. Phys. Commun. 231 (2018) 94.
[4] G. M. Staebler et al., Phys. Plasmas 12 (2005) 102508; ibid. 14 (2007) 055909.
[5] M. Honda and E. Narita, Phys. Plasmas 26, 102307 (2019).
[6] N. Hayashi et al., Proc. Int. Conf. on Fusion Energy (Ahmadabad, India, 2018), PPC/1-1 (2018).
[7] M. Honda and A. Fukuyama, Nucl. Fusion 46 (2006) 580.
[8] P. B. Snyder et al., Phys. Plasmas 16, 056118 (2009).

Reduced models of turbulent transport including zonal ﬂows andtrapped electrons
S. Toda, M. Nunami, and H. Sugama
National Institute for Fusion Science, Oroshi-cho 322-6, Toki, Gifu, 509-5292, JapanThe quantitative prediction of the turbulent transport in toroidal plasmas is one of the mostcritical issues to be solved for realization of fusion energy. Recently, a large number of gyrokinetic
simulations of the turbulent transport in toroidal plasmas have been performed. The gyrokinetic
simulation results for tokamak and helical plasmas have been compared with the experimental
observation results [1]. Since it is known that zonal ﬂows can regulate the turbulent transport,
numerous studies have been done to investigate the eﬃciency of zonal ﬂows to improve plasma
conﬁnement in toroidal devices. In these studies, nonlinear gyrokinetic simulations have been
performed to accurately determine the relation between the turbulent transport level and the zonal
ﬂow amplitude. However, such nonlinear simulations require a huge computational cost for theparameter scan in wide ranges of magnetic ﬁeld conﬁgurations and plasma equilibrium proﬁles. To
reduce the computational cost, the reduced models are proposed. These models can quickly predict
the nonlinear simulation results of the turbulent transport coeﬃcients and ﬂuxes from the linear
simulation results of the instabilities’ growth rates and the zonal ﬂow responses in helical plasmas
under the conditions of adiabatic electrons [2] and kinetic electrons [3]. These reduced models arepresented for the plasmas in the Large Helical Device (LHD), where the ion temperature gradient(ITG) modes are unstable. The turbulent and the zonal ﬂow potential ﬂuctuations can also be
estimated by the reduced model derived from the linear simulation results. To evaluate turbulent
electron particle and heat transport ﬂuxes as well as for treating electromagnetic eﬀects whichbecome important for high-β plasmas, kinetic electrons need to be treated in linear and nonlinear
gyrokinetic simulations. In particular, in helical plasmas, trapped electrons show complicateddrift motions and it is a serious challenge to quantitatively clarify how they impact instabilities,
zonal ﬂows, and turbulent transport. In this study, eﬀects of trapped electrons on zonal ﬂows
and turbulent transport in the LHD conﬁguration are studied. The residue zonal ﬂow level forthe case of kinetic electrons is compared with that in the adiabatic electron condition by thelinear gyrokinetic simulation for the plasmas in the LHD, where the ITG mode is unstable. The
residue level of zonal ﬂows has been studied in tokamak and helical plasmas by the gyrokinetic
simulations. In addition, the zonal ﬂow decay time for the kinetic electron case is compared withthat in the adiabatic electron condition. Next, evaluating the mixing length estimate and the zonal
ﬂow decay time by the linear gyrokinetic simulations, the saturation levels of turbulence and zonal
ﬂows are predicted from the reduced transport models with setting diﬀerent ﬁeld conﬁgurations
and plasma equilibrium proﬁles for which eﬀects of electrons on zonal ﬂows and plasma conﬁnement
are investigated.
The turbulence driven by the microinstabilities and zonal ﬂows in the LHD plasmas are studied,using the gyrokinetic local ﬂux tube code, GKV [4]. In this study, the diﬀusivity and quasilinear
ﬂux models [2, 3] for the ion heat transport are used for the adiabatic electron condition and the
kinetic electron condition. These reduced models are constructed, where the ITG mode is unstable.
The models for the electrostatic turbulent and zonal ﬂow potential ﬂuctuations are also used. The
electron density perturbation δne is given in terms of the electrostatic potential ﬂuctuation byDEi( hδne,k̃x ,k̃ye ϕ̃k̃x ,k̃y − ϕ̃k̃x ,k̃y /Te if k̃y = 0
(1)
=
n0
/Te
if k̃y ̸= 0,eϕ̃k̃x ,k̃ywhere e is the elementary charge, n0 is the background density, Te is the electron temperature,ϕ̃ = ϕ/((Ti /e)(ρi /R)), R is the major radius and ρi (= mi vti /(eB)) isp
the thermal ion gyroradius.
Here, mi is the ion mass, B is the magnetic ﬁeld strength, vti (= Ti /mi ) is the thermal ion
velocity and Ti is the ion temperature. The average along the ﬁeld line is denoted by ⟨· · · ⟩ andk̃x (= kx ρi ) and k̃y (= ky ρi ) are the normalized radial and poloidal wavenumbers, respectively. Only
pure ITG modes are treated by the adiabatic electron response for k̃y ̸= 0. For k̃y = 0, the so-called
modiﬁed adiabatic electron response
is assumed.
Themodels
and 
zonal
 turbulent  for the squared22PP/2 areﬂow potential ﬂuctuations, T = k̃x ,k̃y ̸=0 ϕ̃k̃x ,k̃y
/2 and Z = k̃x ϕ̃k̃x ,k̃y =0
also represented by the linear simulation results, to reproduce the nonlinear simulation results.
The major radii of the plasmas are given by R = 3.75m for the standard ﬁeld conﬁguration and
R = 3.6m for the inward-shifted ﬁeld conﬁguration. The abbreviations SD and IW used in thisstudy stand for the standard and inward shifted ﬁeld
D conﬁgurations,E D respectively.E The linear zonal
ﬂow response function is deﬁned by Rk̃x (t) ≡ ϕ̃k̃x ,k̃y =0 (t) / ϕ̃k̃x ,k̃y =0 (t = 0) . Note that thezonal ﬂow response function for k̃x = 0.25 is used to evaluate the representative values of thezonal ﬂow decay time, because there are peaks of the wavenumber spectra around k̃x = 0.25 in the
nonlinear simulation results. To study the correlation between Rk̃x (t) and the ﬂuctuation of zonal
ﬂows Z̄, the zonal ﬂow decay time is employed. The normalized zonal ﬂow decay time is deﬁned
as τ̃ZF = τZF /(R/vti ). The helical magnetic structure in the inward-shifted ﬁeld conﬁguration
enhances the zonal ﬂow generation.
The linear simulation results related with zonal ﬂows for the modiﬁed adiabatic electron (MAE)
condition of (1) are compared with those for the case of the kinetic electron (KE) response in
the SD and IW. The decay time of zonal ﬂows depends on the magnetic ﬁeld conﬁguration andis independent of the density and temperature proﬁles. When the ﬁeld conﬁguration model wasused, it was analytically shown that the linear response function decays faster due to the presenceof the trapped particles than for the MAE condition. It was found that the residue level for zonal
ﬂows complicatedly depends on the helical ripple and the ratio of the neoclassical polarization due
to toroidally trapped particles to the classical polarization. The linear zonal ﬂow response function
for the case of the KE response is conﬁrmed to decay faster than that for the MAE response in the
SD and IW. The residue level of zonal ﬂows for the case of the KE response is smaller than thatfor the MAE response in the SD at ρ = 0.65. This tendency is the same in the SD at ρ = 0.80.
However, it is clariﬁed that the residue levels are found to be almost the same for the MAE andKE conditions at ρ = 0.65 in the IW. For comparing the zonal ﬂow decay time for the case of
the KE response with that for the MAE condition, the value of τf is set to be 30R/vti for the
MAE and KE conditions. The zonal ﬂow decay time is larger for the MAE condition than for the
KE condition in the SD. However, the zonal ﬂow decay time for the KE condition is found to becomparable to or slightly larger than that for the MAE condition in the IW. The diﬀerence of themagnitude relationship of the residue level for zonal ﬂows or the zonal ﬂow decay time in the SD
and IW is due to the complicated dependence on the helical ripple and the eﬀect of the toroidally
trapped particles.
The nonlinear simulation results, such as the ion heat diﬀusivity, χi as well as the turbulent
and the zonal ﬂow potential ﬂuctuations, T̄ and Z̄, are predicted by the models. The model
predictions for the case of the KE response are compared with those for the adiabatic electron(AE) condition. The values of T̄ for the case of the KE response are found to be larger thanthose for the AE condition in both ﬁeld conﬁgurations, because the trapped electron enhances thegrowth rate of the ITG. Turbulence potential ﬂuctuations in the IW are close to that in the SD.

The models can predict the nonlinear simulation results for T̄ , except for the region ρ < 0.6 in
the IW. The zonal ﬂow potential ﬂuctuation, Z̄ in the IW is found to become larger than thatin the SD, especially for the case of the KE response. The simulation result conﬁrms that zonal
ﬂow potential ﬂuctuation in the IW is larger than that in the SD. The model for Z̄ reproduces
the nonlinear simulation results in the SD. The model does not predict the nonlinear simulation
results especially for the case of the KE response in the IW, because the plasma proﬁle in the SD
is not used for the IW simulation for constructing the model for Z̄. In the SD, the zonal ﬂow
eﬀect is stronger for the the MAE condition than for the KE condition. On the other hand, in theIW, the zonal eﬀect for the case of the KE response is comparable to that for the MAE condition.
The linear simulation results are qualitatively reﬂected in these nonlinear simulation results forthe zonal eﬀect. It is clariﬁed that the reduced models reproduce the nonlinear simulation results.
Even if the linear and nonlinear simulation results in the IW using the plasma proﬁles in the
SD were not adapted for constructing the reduced models for the case of the KE response, thereduced model can reproduce the nonlinear simulation results in the IW for the case of the KE
response. It is found that the transport predicted by the reduced models in the IW is smaller
than that in the SD. Eﬀects of trapped electrons on zonal ﬂows and turbulent transport in two
kinds of LHD ﬁeld conﬁgurations are studied for the representative plasmas, where the ITG mode
is unstable. For the linear simulation results, the residue level and the decay time of zonal ﬂows
for the case of the kinetic electron response is smaller than for the modiﬁed adiabatic electron
response in the SD ﬁeld conﬁguration. However, the residue level and the decay time of zonal
ﬂows for the case of the kinetic electron response are found to be comparable to those for themodiﬁed adiabatic electron response in the IW shifted ﬁeld conﬁguration. For the predictionsof the reduced models, the zonal ﬂow eﬀect on the transport for the case of the kinetic electron
response is weaker than that for the adiabatic electron response in the SD ﬁeld conﬁguration. In
the IW shifted conﬁguration, the zonal ﬂow eﬀect for the case of the kinetic electron response isfound to be comparable to that for the adiabatic electron response. The zonal ﬂow eﬀect on thetransport depends on not only the zonal ﬂow decay time but also the mixing length estimate in
the reduced models. The ion heat diﬀusivities predicted by the reduced model in the IW shifted
ﬁeld conﬁguration are smaller than those in the SD ﬁeld conﬁguration. It is found that the linear
simulation result, such as the zonal ﬂow decay time, is the possible qualitative criterion for zonal
ﬂow eﬀect on the transport. The nonlinear simulation results are well reproduced by the reduced
models for the heat diﬀusivity and the quasilinear ﬂux model, even in the parameter region for thecase of the kinetic electron response in the IW shifted ﬁeld conﬁguration, where the data for theheat transport are not used for constructing the reduced models. This work was partly supportedby the NIFS Collaboration Research Programs, NIFS18KNST129 and NIFS18KNXN363 (Plasma
Simulator), NIFS18KNTT045, and the Collaborative Research Program of Research Institute forApplied Mechanics, Kyushu University, 2019FP-4.
[1] S. Toda et al., Plasma and Fusion Research 14, 3403061 (2019).
[2] M. Nunami, T. -H. Watanabe and H. Sugama, Phys. Plasmas 20, 092307 (2013).
[3] S. Toda et al., Phys. Plasmas 26, 012510 (2019).
[4] T. -H. Watanabe and H. Sugama, Nucl. Fusion 46, 24 (2006).

Flow helicity of wavy plasma turbulence
Y. Kosuga[1,2][1] Research Institute for Applied Mechanics, Kyushu University, Fukuoka, Japan[2] Research Center for Plasma Turbulence, Kyushu University, Fukuoka, Japan
Helicity of flow field, defined as v · ω where v is the velocity and ω = ∇ × v is thevorticity, is an important quantity to understand the dynamics of turbulent plasmas.
For example, turbulence with finite flow helicity can transport momentum to generatelarge scale flows. Turbulence with finite helicity also can generate large scale magnetic
field, which can be modelled via the famous alpha term in the mean field theory of themagnetic field evolution. Moreover the generation of large scale helical flows is reported
from theory[1] and numerical simulation[2]. Flow helicity plays an important role toquantify the resultant flow structures.
In this work[3], we have analyzed how the helicity of turbulent plasmas is determined,by using drift wave turbulence and parallel velocity gradient driven (PVG) turbulence.
First of all, we identify relevant symmetry breaking to have the flow helicity to be finite.
Generally, the helicity of the flow field is finite when the reflectional symmetry is broken in
the system. In the case of magnetized plasmas, we find that the breaking of the reflectional
symmetry is closely related to the breaking of the symmetry along the magnetic field. In
this regard, PVG turbulence is a typical example of the turbulence with finite flow helicity,
since the growth is asymmetric with regard to the direction of the magnetic field. For the
given parallel flow shear and the poloidal mode number, fluctuation with the particularsign of the parallel wave number is excited. This translates into the helicity of PVG
turbulence asccs /Ln ∑ 2 2 eϕ̃k⟨ṽ · ω̃⟩ = 2 ⟨ṽz ∇2⊥ ϕ̃⟩ = 2cs ωcik ρB⟨vz ⟩′ k ⊥ s Te

2

(1)

Thus the helicity of PVG turbulence, especially the entire sign, is determined by ∝⟨n⟩′ ⟨vz ⟩′ B. Then flipping the magnetic field or the parallel flow direction, the sign ofthe helicity flips. Pictorially, turbulence may be filled with the lefthand sided helix, whileby flipping these sings the turbulence can be filled with the righthanded helix.

Then the dynamical evolution of fluctuation helicity is formulated. Especially weemphasized the link of the evolution of fluctuation helicity to the transport fluxes. For
PVG turbulence, the relevant flux to drive the fluctuation helicity evolution is tied tothe momentum flux of parallel flows. Note that the parallel momentum flux requires thesymmetry breaking to have a finite value. The finiteness of this flux is tied to the finiteness of the fluctuation helicity, which is a plausible result. Then the obtained evolution
equation is used to diagnose the simulation data obtained from the NLD code. A steadystate balance is achieved by balancing the damping of helicity and the parallel momentum
flux. The resultant expression has the similar dependence on the sign as Eq.(1).
Stimulating discussion with P.H. Diamond, K. Itoh, S. Inagaki, F. Kin, is acknowledged. This work was supported by Grants-in-Aid for Scientific Research of JSPF ofJapan (JP15H02155, JP17H06089, JP18K03578), the Joint research project in RIAM.

References[1] Y. Kosuga, S.-I. Itoh and K. Itoh, Phys. Plasmas 24 032304 (2017)[2] M. Sasaki, N. Kasuya, K. Itoh, S. Toda, T. Yamada, Y. Kosuga, Y. Nagashima, T.
Kobayashi, H. Arakawa, K. Yamasaki, A. Fujisawa, S. Inagaki, and S.-I. Itoh, Phys.
Plasmas 24 112103 (2017)[3] Y. Kosuga, M. Sasaki, Z.B. Guo Phys. Plasmas 27 022303 (2020)

Eﬀects of Electron Cyclotron Heating on the Toroidal Flow in
HSX Plasmas
Y. Yamamoto, S. Murakami, S.T.A. Kumar∗ , J.N. Talmadge∗ , K.M. Likin∗ , D.T. Anderson∗ ,Department of Nuclear Engineering, Kyoto Univ., Nishikyo, Kyoto 615-8540, Japan∗ Department of Electrical and Computer Engineering, University of Wisconsin-Madison,
Madison, WI 53706, USAe-mail: yamamoto y@p-grp.nucleng.kyoto-u.ac.jp
The toroidal flow related to electron cyclotron heating (ECH) is investigated in HSX plasmas. Radialdiﬀusion of energetic electrons by ECH produces a canceling return current, which then generates
a jr × B torque that can play an important role in the toroidal rotation in ECH plasma. We
investigate the energetic electron distribution by ECH by applying GNET code, which can solve the5D drift kinetic equation for the energetic electron. We evaluate the jr ×B torque and the collisional
torque due to the friction of the toroidal drift motion of the energetic electrons. As a result, we
obtained a significant torque due to the ECH and found that a larger torque is obtained in the
Mirror configuration than in the quasi-helically symmetric (QHS) configuration in HSX. Solving the
momentum balance equations and Ampère’s law with the jr ×B torque, we evaluate the toroidal flow
velocity and compare simulations with HSX experiments. The obtained flows have good agreementwith experimental ones.
Introduction
Recently, spontaneous toroidal flows have been observed inelectron cyclotron heating (ECH) plasmas in many tokamak andhelical devices such as JT-60U, LHD and HSX. To clarify theunderlying mechanism, many experimental [1] and theoreticalstudies have been undertaken. The eﬀects of the magnetic configuration on plasma flow are intensively investigated in HSX,
where two typical magnetic configurations are considered. One
is the Quasi-Helically Symmetric (QHS) configuration, which
has a quasi-helical symmetry in |B| and is dominated by the[m, n] = [1, 4] mode. The other is the Mirror configuration, wherea set of auxiliary coils adds toroidal mirror terms, the (0, 4) and(0, 8) modes, to the magnetic field spectrum to break the helicalsymmetry. The QHS configuration has a parallel neoclassical viscosity that is smaller than that of the Mirror configuration, so we Figure 1: The total helical forces
expected that the toroidal flow velocity in the QHS configuration including j × B and collisionalrwould be larger than that of the Mirror configuration. However torques in three configurations.
smaller toroidal flow was observed in the QHS configuration. It
has not been understood well yet.
Simulation model
ECH can drive the radial electron current je due to the radial motion of suprathermal electrons[2]. The net current in the steady state should be canceled to maintain the quasi-neutrality, so thereturn current, jr (= −je ) , must flow by the bulk ions by ambipolar condition. Therefore, thebulk plasma feels the jr × B torque due to the return current. On the other hand, the suprathermal
electrons drift toroidally due to the precession motion. During the slowing down of the suprathermal
electrons, they transfer their obtained momentum to the bulk plasma due to collisions.
In this study, we investigate the behaviors of energetic electrons by ECH, which can generate
the radial current making the jr × B torque in the HSX plasmas. Also, we evaluate the collisional
torques, by the collision between energetic electrons and bulk plasma. We apply the GNET code,1which can solve a linearized drift kinetic equation for energetic electrons by ECH in 5-D phasespace[2]
∂δf
∂δf
∂δf
+ (⃗vd + ⃗v∥ ) ·+ v̇ ·− C(δf ) − L(δf ) = S ql (fMax ).
∂t
∂⃗r
∂⃗v
Solving the momentum balance equations and Ampère’s law withthe radial electron flux and torques, we evaluate the toroidal flow
velocity and compare them with experimental one.
Results
In the perfectly symmetric configuration, the forces in thesymmetry direction cancel each other. As seen in Fig.1, the jr ×B
torque and the collisional torque cancel each other and the component parallel to the helically symmetry direction is very small
in the perfectly helically symmetric configuration[3]. However
non-symmetric magnetic modes enhance the radial current je .
Thus, even in the QHS configuration the jr × B torque is dominant and there is a net force in the symmetry direction due toother small non-symmetric modes. The force in Mirror configuration is more than twice as large as that in QHS configuration withthe same input power. The collisional torque is so small that we
can ignore it in QHS and Mirror configurations. Solving the momentum balance equations with jr × B torque, the obtained flow
velocity is shown in Fig.2. Here the absorption power calculatedby ray-tracing code is 24kW in QHS configuration and 16kW in
Mirror configuration. The obtained flow in QHS configuration isstrongly peaked, but the toroidal flow in Mirror configuration islarger than that in QHS configuration totally.

Figure 2: The obtained toroidal
flow in QHS and Mirror configurations. The lines are simulation results and the dots are experimental ones.

[1] S.T.A. Kumar, et al., Nucl. Fusion 60 (2018) 054012.
[2] S. Murakami, et al., Nucl. Fusion 40 (2000) 693.
[3] Y. Yamamoto, et al., Plasma and Fusion Research 14 (2019) 3403105[4] M. N. Rosenbluth and F. L. Hinton, Nucl. Fusion 36 (1996) 55

2

Impacts of plasma rotation and ion diamagnetic drift on MHD
stability in a QH-mode plasma

N. Aiba1*, X. Chen2, T. H. Osborne2, S. R. Haskey3, K. H. Burrell2, P. B. Snyder2
1 National Institutes for Quantum and Radiological Science and Technology, Naka, Ibaraki
311-0105, Japan
2 General Atomics, PO Box 85608, San Diego, CA 92186-5608, USA
3 Princeton Plasma Physics Laboratory, PO Box 451, Princeton, NJ 08543-0451, USA* Email address of submitting author: aiba.nobuyuki@qst.go.jp
In H-mode regime in tokamaks, edge localized modes (ELMs) often appear and induce large
heat load to divertors. Since the heat load is unacceptable for future large reactors like ITER
and DEMO, it is necessary to suppress and/or mitigate the ELMs. Quiescent H-mode (QHmode) is one of the promising candidates realizing ELM suppression and high confinement
performance with ITER and DEMO relevant plasma parameters [1]. One of the importantcharacteristics of the QH-mode is that edge harmonics oscillations (EHOs) are observed,
though ELMs disappear. Since QH-mode can be obtained when bootstrap current and rotation
shear near plasma surface are large, the EHOs have been recognized as a rotation shear
destabilized current-driven MHD (peeling) mode. The hypothesis has been validated by linear
analyses with the ideal [2] and resistive [3] MHD models.
In addition, the ion diamagnetic drift frequency, ω*i, can be comparable to the plasma rotation
frequency[3], hence, the 𝐸𝐸 × 𝐵𝐵 rotation has been recognized as a strong candidate responsible
for entering QH-mode, the rotation which is composed of the sum of plasma fluid rotation andthe ion diamagnetic drift [4]. The MHD analysis with the 𝐸𝐸 × 𝐵𝐵 rotation as already performedin [3], however, in this case, the ion diamagnetic drift was treated as the single fluid MHD
rotation. The stabilizing effect by the ion diamagnetic drift on pressure-driven ballooning
mode in a static plasma has long been recognized as discussed in [5], and the effect has beenobtained by treating the fluid rotation and the ion diamagnetic drift independently.
In this study, we analyze numerically the MHD stability in DIII-D QH-mode plasmas byconsidering plasma rotation and ω*i effects simultaneously with the linear extended MHD
stability code MINERVA-DI [6]. This code solves the extended Frieman-Rotenberg (F-R)
equation∂𝛏𝛏𝜕𝜕𝝃𝝃⊥𝜕𝜕 2 𝝃𝝃
𝜌𝜌0 2 + 2𝜌𝜌0 �𝑽𝑽0,𝑀𝑀𝑀𝑀𝑀𝑀 ⋅ ∇� + 𝜌𝜌0 �𝑽𝑽0,∗𝑖𝑖 ⋅ ∇�= 𝐹𝐹(𝝃𝝃),
𝜕𝜕𝜕𝜕
∂t
𝜕𝜕𝑡𝑡
∇ ⋅ 𝝃𝝃=0
with the definition of velocity vectors as
𝑉𝑉∥ 𝑩𝑩
𝑽𝑽𝑀𝑀𝑀𝑀𝑀𝑀 = 𝑽𝑽⊥ +,|𝐵𝐵|
1
𝑽𝑽∗𝑖𝑖 =
𝑩𝑩 × ∇𝑝𝑝𝑖𝑖 ,̅
𝑒𝑒𝑍𝑍𝑁𝑁|𝑩𝑩|2
where 𝐹𝐹(𝜉𝜉) is the linearized force operator derived from the diamagnetic MHD model [6].
Here 𝝃𝝃 = 𝝃𝝃⊥ + 𝜉𝜉∥ 𝑩𝑩/|𝐵𝐵| is the Lagrangian displacement vector defined with the linearized
Eulerian velocity, 𝜌𝜌 is the plasma mass density, 𝑒𝑒 is the quantum of electricity, 𝑍𝑍̅ is the ion
mean charge satisfying 𝑁𝑁𝑒𝑒 = 𝑍𝑍̅𝑁𝑁𝑖𝑖 , 𝑁𝑁𝑒𝑒 (𝑁𝑁𝑖𝑖 ) is the electron (ion) number density, 𝑝𝑝𝑖𝑖 is the ionpressure, and the subscript 0 indicates the equilibrium quantity.
The analyzed DIII-D plasma is #153440, which has the EHO with the large n=2 component inmagnetic fluctuation, where n is the toroidal mode number. The parameters of this plasma are
the toroidal magnetic field |𝐵𝐵𝑡𝑡 |~1.91𝑇𝑇, the plasma current �𝐼𝐼𝑝𝑝 �~1.5𝑀𝑀𝑀𝑀, 𝐵𝐵𝑡𝑡 and 𝐼𝐼𝑝𝑝 areclockwise from the top view of the torus, the ellipticity 𝜅𝜅~1.85, and the triangularity 𝛿𝛿~0.5.
As discussed in [3,4], the rotation shear in not the toroidal rotation of impurity (carbon) but
the toroidal rotation associated with the 𝐸𝐸 × 𝐵𝐵 drift, ω𝐸𝐸×𝐵𝐵 , is thought to be responsible forentering the QH-mode, hence, the MHD stability is analyzed by taking into account Ω𝐸𝐸×𝐵𝐵 .
Note that such the analysis can be realized by setting 𝑽𝑽0,𝑀𝑀𝑀𝑀𝑀𝑀 = 𝑅𝑅ω𝜙𝜙 − 𝑟𝑟ω𝜃𝜃 𝐵𝐵𝑡𝑡 /|𝐵𝐵𝜃𝜃 | =𝑅𝑅ω𝑣𝑣×𝐵𝐵 , because the ion diamagnetic drift frequency in ω𝐸𝐸×𝐵𝐵 is explicitly included as 𝑽𝑽0,∗𝑖𝑖 inthe extended F-R equation. Here Ω𝜙𝜙 (Ω𝜃𝜃 ) is the rotation frequency in the toroidal (poloidal)
direction, 𝑟𝑟 is the plasma minor radius, and 𝐵𝐵𝑝𝑝 is the poloidal magnetic field. The ω𝜙𝜙 and 𝜔𝜔𝜃𝜃
profiles are determined by using the measured rotation profiles of carbon, ω𝜙𝜙,𝐶𝐶 and ω𝜃𝜃,𝐶𝐶 , andthe corresponding ω𝑣𝑣×𝐵𝐵 is named as ω𝑣𝑣×𝐵𝐵,𝐶𝐶 . The ion diamagnetic drift frequency 𝜔𝜔∗𝑖𝑖 profile
is evaluated with both the deuterium and carbon pressure profiles with the assumption that thetemperature is the same as each other. The equilibrium profiles of the plasma are shown in
Fig.1.

Fig.1 Profiles of the DIII-D 153440/1.724sec. plasma; (left) 𝑇𝑇𝑖𝑖 , 𝑇𝑇𝑒𝑒 , 𝑁𝑁𝑖𝑖 , 𝑁𝑁𝑒𝑒 ,(center) 𝑃𝑃0 , 𝑞𝑞, (right) 𝜔𝜔𝜙𝜙,𝐶𝐶 , 𝜔𝜔𝜃𝜃.𝐶𝐶 , 𝜔𝜔𝑣𝑣×𝐵𝐵,𝐶𝐶 , 𝜔𝜔∗𝑖𝑖

Figure 2 shows the stability diagram of the plasma on the 𝑗𝑗𝑝𝑝𝑝𝑝𝑝𝑝 − 𝛼𝛼 plane, where 𝑗𝑗𝑝𝑝𝑝𝑝𝑝𝑝 is theflux averaged current density inside edge pedestal, and is the flux averaged current density
inside edge pedestal, and 𝛼𝛼 is the normalized pressure gradient. When analyzing the stability
without the 𝜔𝜔∗𝑖𝑖 effect, the stability boundary indicates that the current driven MHD mode
(kink/peeling mode) is destabilized by rotation because the mode can be unstable with small
amount of 𝑗𝑗𝑝𝑝𝑝𝑝𝑝𝑝 as shown in Fig.2(a). The trend is consistent with the previous results
discussing the impact of sheared rotation on kink/peeling mode stability. However, when
considering the 𝜔𝜔∗𝑖𝑖 effect simultaneously, the kink/peeling mode in rotation plasmas requires
large amount of 𝑗𝑗𝑝𝑝𝑝𝑝𝑝𝑝 for destabilization as shown in Fig.2(b). Namely, the rotation stabilizes

Fig.2 MHD stability diagram of the DIII-D 153440/1.725sec. QH-mode
plasma on the 𝑗𝑗𝑝𝑝𝑝𝑝𝑝𝑝 − 𝛼𝛼 diagram; the target expresses the equilibrium values
observed in experiment, and the numbers on the thagram show the toroidal
mode number of the MHD mode determining the stability boundary; (left)
without 𝜔𝜔∗𝑖𝑖 , (right) with 𝜔𝜔∗𝑖𝑖 .

the king/peeling mode. The trend is clearly observed when varying the rotation profile bymultiplying coefficient C to 𝜔𝜔𝑣𝑣×𝐵𝐵,𝐶𝐶 ; the result is shown in Fig.3. As shown in this figure, thestability boundary of the peeling modes slightly moves downward by increasing 𝜔𝜔𝑣𝑣×𝐵𝐵,𝐶𝐶 ,because the mode is destabilized. On the other hand, when including both rotation and the 𝜔𝜔∗𝑖𝑖
effects, the boundary moves upward; namely the mode is stabilized by rotation.
These results imply that the rotation effect on kink/peeling mode stability depends on theexistence of the ion diamagnetic drift effect. Since the ion diamagnetic drift at edge pedestal
region in QH-mode plasmas is usually large, both effects could play key roles for entering
QH-mode. To validate this numerically, nonlinear simulation including both the effects isnecessary to obtain QH-mode plasma numerically as a steady state with EHOs. This remains
as our future work.

Fig.3 Impact of rotation frequency, 𝐶𝐶𝜔𝜔𝑣𝑣×𝐵𝐵,𝐶𝐶 , on the kink/peeling modestability boundary on the 𝑗𝑗𝑝𝑝𝑝𝑝𝑝𝑝 − 𝛼𝛼 diagram; (left) without 𝜔𝜔∗𝑖𝑖 , (right)
with 𝜔𝜔∗𝑖𝑖 .

This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy
Sciences, using the DIII-D National Fusion Facility, a DOE Office of Science user facility, under Awards DE-FC02-04ER54698.
DIII-D data shown in this paper can be obtained in digital format by following the links
at https://fusion.gat.com/global/D3D_DMP.

[1] K. H. Burrell et al., PRL 102, 155003 (2009).
[2] T. H. Osborne et al., J. Phys. Conf. Ser. 120 012014 (2008).
[3] X. Chen te al., Nucl. Fusion 56 076011 (2016).
[4] A. Garofalo et al., Nucl. Fusion 51 083018 (2011).
[5] M. S. Chu et al., Phy. Rev. Lett. 41 247 (1978) etc..
[6] N. Aiba, PPCF 58, 045020 (2016).

Nonlinear dynamics of resistive ballooning instability
in PLATO tokamak plasma
S. Tomimatsu , N. Kasuya , M. Sato , A. Fukuyama , M. Yagi , Y. Nagashima ,
1

1,2

3

and A. Fujisawa
1

4

5

1,2

1,2

IGES, Kyushu University, RIAM, Kyushu University, NIFS, Kyoto University, QST
2

3

4

5

1. Introduction
In order to realize a fusion reactor, it is necessary to sustain stable confinement of high-temperatureplasmas. However, in a torus plasma, there are various instabilities that have their characteristic spatialtemporal scales (e.g. MagnetoHydroDynamics (MHD), ion and electron scales), and these instabilities
make it difficult to confine the plasma [1,2]. MHD instabilities due to the pressure gradient and the plasma
current cause plasma collapse phenomena [3-5]. Therefore, understanding the characteristics of MHD
instabilities is essential for stable control of plasma performance.
In order to understand the physical mechanism of plasma dynamical phenomena, laboratory size
tokamak PLATO is now under construction in Kyushu university to make detailed measurements of plasma
perturbations by using various experimental diagnostics [6]. The flexible design of PLATO tokamak makes
it possible to observe dynamics of 3-D turbulence structures. The expected plasma parameters of PLATO
tokamak parameters are following; major radius 𝑅 = 0.7m, minor radius 𝑎 = 0.2m, plasma volume 𝑉* =19 3
0.9m3, electron temperature 𝑇- = 0.2keV, density 𝑛 = 1.0×10 m .
In this study, the target is MHD instabilities with macro scales, and we carry out MHD simulations byusing PLATO tokamak parameters to represent dynamics of MHD instabilities in a tokamak plasma.
Ballooning and kink modes are induced depending on the safety factor q profile and the steep gradientregion, and pressure flattening processes accompanied with competition of instabilities are investigated in
the nonlinear dynamics.

2. Equilibrium in PLATO
The plasma equilibria considering the experimental conditions areevaluated by using the module EQU in integrated transport code TASK [7].
TASK/EQU calculates equilibria with free boundary conditions in existenceof external vertical coils [8]. Figure 1 shows the magnetic configuration
obtained by using the set of external coil currents. In this calculation, only D,
H and Q coils are used for control parameters. This calculation gives toroidal
magnetic field 𝐵2 = 0.3T, plasma current 𝐼* = 40kA, central-beta 𝛽2 = 2%,
central-q 𝑞<=>? = 0.7 , surface-q 𝑞@ABC = 9.8 , ellipticity 𝜅 = 1.6
triangularity 𝛿 = 0.3.

3. MHD simulation

Q
F2
H
F1
C
D

0.8

0.4


Z
0.0
-0.4

-0.8


0.0 
0.4 
0.8 1.2

 1.6

R

Fig.1Magnetic configuration obtained by TASK/EQU.
The configuration has up-down symmetry.

For simulations of MHD instabilities in PLATO tokamak, MIPS code [9] is used. MIPS code solves
the MHD equations described below in the cylindrical coordinates (𝑅, 𝑍, 𝜙);

𝜌

𝜕𝜐
𝜐T
4
= ρ𝜔×υ − ρ∇− ∇𝑝 + 𝚥×𝐵 + ∇ 𝜈𝜌 ∇ ∙ 𝜐
𝜕𝑡
2
3

𝜕𝜌
= −∇ ∙ 𝜌𝜐 ,𝜕𝑡− ∇× 𝜈𝜌𝜔 ,

(1)

(2)

\]

= −∇×E,(3)

𝐸 = −υ×B + 𝜂𝚥,

(4)

1
∇×𝐵,
𝜇c

(5)

\2

𝚥=

𝜔 = ∇×𝜐.

(6)

𝜕𝑝= −∇ ∙ 𝑝𝜐 − 𝛾 − 1 𝑝∇ ∙ 𝜐 + 𝜒∇T 𝑝 − 𝑝-g𝜕𝑡4+ 𝛾 − 1 𝜈𝜌𝜔 T + 𝜈𝜌 ∇ ∙ 𝜐 T + 𝜂𝚥 ∙ 𝚥 − 𝚥-g ,3(7)Here, 𝜌, p, 𝜐, 𝐵, 𝐸, 𝚥 are the mass density, plasma pressure, velocity, magnetic field, electric field andcurrent density, 𝑝-g , 𝚥-g are the equilibrium pressure and current, respectively. 𝜈, 𝜂, 𝛾, 𝜒 are viscosity,
resistibility, adiabatic constant and thermal diffusion coefficient, and 𝜈 = 𝜂/𝜇c = 10ij 𝑣l 𝑅c , 𝜒 =10im 𝑣l 𝑅c are used for the following calculations, where 𝑣l is the Alfvén velocity, 𝑅c = 0.7 m is the
typical major radius. The time t is normalized by the Alfvén time.
The MHD equilibria used in the simulations are constructed by TASK/EQ module which calculates
equilibria with a fixed boundary condition. The magnetic surface of q = 5.5 of the equilibrium calculatedby the TASK/EQU module shown in Fig.1 is used for the boundary condition of the TAKS/EQ module.
Figure 2 shows the time evolution of the pressure perturbation energy for qaxis = 0.8, where the pressure
is normalized by B02/µ0. The total pressure perturbation energy, which is integrated for the whole plasma
volume, increases monotonically in the linear phase, and show nonlinear saturation after the phase. Its
Fourier decomposition to the toroidal components is also shown in Fig. 2. The initial growth rate of the n= 1 component is smaller than that of the n = 14 component, where n is the toroidal mode number. The
growth rate of the n = 1 component increases at t = 300, and at t = 500 its magnitude becomes larger thanthose of the other modes. The pressure components in the poloidal cross-section are shown in Figs. 3. In
Fig. 3 (t = 300), a ballooning mode structure whose perturbation component exists in the low field side ofplasma [10,11] can be seen with finite poloidal mode numbers. Next, we conducted comparison of the mode
spectra in the linear phase (t = 200) and nonlinear phase (t = 400). In the linear phase, the dominant modes
have the poloidal mode numbers of m = 19 around 𝜌 = 0.65, where the pressure gradient is larger. Herem is the poloidal mode number. On the other hand, in the nonlinear phase, dominant modes exist around𝜌 = 0.45. This is because the pressure profile near 𝜌 = 0.65 is flattened by the instabilities, and thesteeper gradient region is shifted to the inner part near 𝜌 = 0.45. In this way, the mode excited region is
expanded in accordance with the flattening process. After the flattening by the ballooning modes, the n = 1
mode induces sudden change of the pressure profile in the central region as in Fig. 3 (t = 500 and t = 600).
A m/n =1/1 kink mode at the q = 1 surface [12] gives the collapse near the center. In this way, both ofballooning and kink modes are excited in this simulation.

Total
n=1
n=14

10-4

×10$%

0.4

10-6
10-8

10-10

0

200

t

400

t=0

t = 300

t = 500

1.2

t = 600

1.0

0.2

0.8

Z 0.0

0.6

-0.2

0.4

-0.4

0.2
0.4 0.6 0.8
R600

Fig. 2 Time evolutions of the pressure perturbationenergy.

1 0.4 0.6 0.8
R

1 0.4 0.6 0.8
R

1 0.4 0.6 0.8
R

10

Fig. 3 Time evolution of the pressure profile in thepoloidal cross-section.

4. Characteristic of the instabilities

Normalized pressure (m,n)=(0,0)

In this section the characteristics of the instabilities are investigated by changing plasma parameters;
plasma beta, resistivity and safety factor profile. First, the dependency of the ballooning instability on theplasma beta is explained. The linear growth rate becomes larger, as the plasma beta is increased from 1%
to 6%. Therefore, it is revealed that this instability is strengthened by the larger pressure gradient. Next, thedependency on the resistivity is examined. We carry out MHD simulations by changing 𝜂 from 10in to
10io . The linear growth rate becomes larger, as 𝜂 is increased in the range 𝜂 > 10ip . Note that 𝜂 = 10ij
is used for the nonlinear simulation in the previous section.
In the MHD simulation in Sec. 3, both the ballooning and kink modes appear. The excited mode canbe controlled by changing the q profile.
Different q profiles are used as the initial
×10$%
(b)
(a)
1.2
conditions. Figures 4 show comparison oft = 175
t = 250t = 175
t = 250t = 275
t = 225
t = 400
the time evolutions of the pressure profilet = 300
0.8
t = 350
with different q profiles. In the Fig.4 (a)t = 400
with qaxis = 1.3, there is no q = 1 surface in
t = 175
0.4
t = 250
the plasma, so the kink mode is not excitedt = 275
t = 400
and the ballooning mode is dominant.
0 0 0.2 0.4
0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
When the ballooning mode is dominant,normalizedradiusnormalized radius
there is no sudden change in the central
Fig.4 Time evolutions of the pressure profile for the (a) ballooningpressure as in Fig. 4(a). On the other hand,dominant case (qaxis = 1.3) and (b) kink dominant case (qaxis = 0.5).
in the Fig.4 (b) with qaxis = 0.5, since qaxis issmaller than the case with qaxis = 0.8, the kink mode is more unstable, which causes sudden change in thecentral pressure as shown in Fig. 4(b). Therefore, the flattening process of the pressure profile is drastic. It
is revealed that the flattening process is changed by dominant instabilities.

5. Conclusion
In this study, the MHD simulations for PLATO plasmas were carried out by the combination of TASK
and MIPS codes. The plasma equilibrium was evaluated with the vertical coil configuration of PLATO by
using the free boundary equilibrium code TASK/EQU. The equilibrium was introduced into the MHD
simulations and nonlinear saturation dynamics was calculated by using MIPS code. The ballooning andkink modes can both become unstable for the PLATO parameters. The parameter dependencies show theircharacteristics. The ballooning modes become unstable at the region with steep pressure gradient, whichgive rather mild flattening of the pressure profile on the low field side of the plasma. On the other hand, thekink modes with m/n = 1/1 become unstable near the center of the plasma, and give drastic pressure collapseof the core.

Normalized pressure

Pressure perturbation energy

10-2

References[1] A. M. Dimits, et al., Phys. Rev. Lett. 77, 1 (1996)[2] R. J. Fonck, et al., Phys. Rev. Lett. 70, 24 (1993)[3] K. L. Wong, et al., Phys. Rev. Lett. 85, 5 (2000)[4] Y. Liang, Phys. Rev. Lett. 98, 265004 (2007)[5] P. B. Synder, et al., Phys. Plasmas 9, 2037 (2002)[6] A. Fujisawa, AIP Conf. Proc. 1993, 020011 (2018)[7] http://bpsi.nucleng.kyoto-u.ac.jp/task/[8] K. Matsuoka, Reports RIAM, Kyushu Univ. 141, 51 (2011)[9] Y. Todo, et al., Plasma Fusion Res. 5, S2062 (2010)[10] J. P. Freidberg, Ideal MHD (Cambridge University Press, 2007)[11] J. Wesson, Tokamaks (Oxford University Press, 2011)[12] R. Y. Dagazian and P. H. Rutherford, Phys. Fluids 16, 1894 (1973)[13] M. Yagi, et al., Plasma Fusion Res. 2, 025 (2007)

Particle transport model in plasma peripheral region
M. Yagi
National Institutes for Quantum and Radiological Science and Technology,
Rokkasho Fusion Institute, Rokkasho, Aomori, Japan
1. Introduction
To understand particle pinch mechanism in the peripheral region of Tokamak is important
issue to supply of D-T fuels in DEMO reactor. The inverted density gradient appears in ahollow density profile when gas-puff or pellets are injected. Using a delta-f electromagnetic
gyrokinetic simulation code dFEFI[1], we investigated a local dynamics of particle transport
in the case with inverted density gradient. It is shown that in the nonlinear phase, the lower
wave number mode becomes dominant, which produces an inward particle flux[2]. With ahigh electron collisionality, electron temperature fluctuations associated with the effects offinite electron thermal conductivity can produce a phase shift between the density andelectrostatic potential fluctuations. The phase shift causes an inward pinch in the ITG-like
mode, so called ‘ion-mixing mode’[3]. On the other hand, we have also observed that thesimilar phenomena in the 5-field reduced MHD simulation with density hump where ITG is
not taken into account[4]. In this work, we investigate the electron drift wave (EDW) withinverted density gradient. It is shown that in this situation, it could drive inward particle flux
rather than the ion-mixing mode.
2. Electron drift wave
The electron response is modeled as1with 𝛿

𝑖𝛿 ,

𝑘 𝜔

(1)
1.5𝜂 𝜔∗ ⁄𝜔 , 𝑘𝜔∗1.76, 𝜔𝑘∥ 𝑣

⁄𝜈 . According to Ref.[5,6],

we obtain
1

∥

𝑖𝛿

𝑏

∗

1

0 .

(3)

Then, the quartic equation is obtained as
𝑖

Ω

1

𝑏

𝑖

where

𝐺

2𝜖 𝐾, 𝜉Ω1𝑘∥ 𝑐 ⁄𝜔∗ , 𝑆

𝑏𝐾

2𝜖 Ω

𝐾𝜉 , 𝜔𝜉𝑘𝜔 ⁄𝜔∗ , 𝜂̂𝐺 Ω𝑆3⁄ 2 𝜂

0,2⁄3

represents the driving term of toroidal ITG mode and 𝑆, the one of slab ITG mode.

(4)
.

𝐺

The particle flux is given by
Γ

𝑛 ∑

If 𝛿

0, Eq.(5) implies the inward particle flux. In Eq.(4), we shall consider the electron

𝑅𝑒 𝛿 .

(5)drift wave with the inverted density gradient, namely, 𝜔∗

0, 𝜉

0, 𝐺

0, Eq.(4) reduces to the quadratic equation as

0, 𝑆

0, 𝜖

Ω

𝑖𝜔

In the limit of𝜔

0. Setting 𝑏𝜔∗𝑖Since 𝜔

𝜂̂ Ω

𝑖𝜔

0 .

𝑖𝜔

𝜂̂

≫ 4𝜔 , the unstable solution is given by

𝜔∗

∗

𝜔∗

(6)

.

(7)

0, the imaginary part is positive. Then,
𝑘

𝑅𝑒 𝛿

∗

0 ,

(8)

which implies the inward particle flux is driven by the electron drift wave with the inverted
density gradient. This model is suitable for semi-collisional regime where the passing electron
plays a role. On the other hand, if the collision frequency enters the banana regime, thetrapped electron may play a similar role.
3. Effect of dissipative trapped electron mode on EDW
We consider the dissipative trapped electron mode (DTEM) on the electron drift wave
assuming 𝜈 ~𝜔 . The electron response is modeled as
1

𝑖𝛿

𝑖𝛿

.

(9)

For the DTEM, the nonadiabatic electron response is given by[7]
𝛿

√

/

√

∗

Ω

1

𝜂

(10)

For Ω~1, it is found that the particle flux is the outward, namely,
Γ

𝑛 ∑

𝑅𝑒 𝛿

√
√

𝑛 ∑

𝜔∗

0 .

(11)

DTEM effect on inverted density gradient driven EDW is easily taken into account of thenonlocal dispersion relation given by Eq.(11) in Ref.[6],
Ω

𝜖∗ Ω

1 1

𝑖𝛿

𝜈∗ Ω Ω

𝜂̂

0.

It implies that DTEM has a stabilizing effect on EDW through the factor 1

(12)
𝛿 .

4. Summary and discussion
We investigate the EDW in the case with the inverted density gradient and normal
temperature gradient in tokamak plasma peripheral region. The local dispersion relation

Eq.(7) shows the electron temperature gradient and the parallel electron thermal conductivity
drives the instability in the case with 𝜔∗

0 which produces the inward particle flux. We

also investigate the effect of DTEM on EDW. It is found that DTEM has a stabilizing effecton EDW through the factor 1

𝛿and weaken the inward particle flux.

To transport the tritium fuel in the core region, we need another mechanism which produces
the inward particle flux. The fluid model which describes ITG/TEM turbulence in thecollisionless regime was proposed by Nordman et al.[8]. In their paper, they investigated theparticle pinch and thermal pinch regime in the phase space of 𝜀 , 𝜂 . There exists theregion where the particle pinch occurs. For this region, ITG is unstable and coupled with
TEM. The growth rate is given by[8],
𝛾

1
1

√𝜖
𝜖 1
√𝜖

𝜂

7
𝜖
4

1

𝜖
4

To couple the plasma core with the edge, we shall consider the phase space
𝜈 ,𝐿 ,𝐿,𝐿

instead. For future work, we could try the Gedanken experiments, namely,(1) assuming quasi-linear flux, we search the region which gives particle pinch and thermal
pinch in parameter space, making phase space map, then, (2) learning steady state profilechanging heat source and particle source to approach such region by the neural network, then,we could make particle control by the neural network.
Acknowledgments
The author would like to thank Drs. N. Miyato, A. Matsuyama, N. Kasuya and T. Takizuka
for useful comments.
References[1] B. Scott, Phys. Plasmas A17 (2010)102306[2] K. Miki, N. Miyato, B. Scott and M. Yagi, 26th International Toki Conference, NIFS,
2017/12[3] B. Coppi and C. Spight, Phys. Rev. Lett. 41 (1978) 551[4] M. Yagi et al., IAEA FEC 2016, Kyoto (2016) THP3-21.
[5] F. Romanelli, Phys. Fluids B 1 (1989) 1018.
[6] M. Yagi, et al., Proceedings of 16th Burning Plasma Simulation Initiative Meeting, Nov.
29-30, 2018, Kyushu University, Kasuga, Japan.
https://www.riam.kyushu-u.ac.jp/sosei/bpsi/2018/BPSI2018_proc.pdf[7] R. Romaneli and S. Briguglio, Phys. Fluids B 2 (1990) 754.
[8] H. Nordman et al., Nucl. Fusion 30 (1990) 983.

Study of multiple impurity seeding effect using SONIC integrated
divertor code for JT-60SA plasma prediction
Shohei YAMOTO1), Kazuo HOSHINO2), Tomohide NAKANO1), Nobuhiko HAYASHI1)
1)

QST, 2)Keio Univ.

1. Introduction
The JT-60SA tokamak is now being constructed and is expected to ignite the first plasma at 2020 andto perform ITER-supporting and -complementing experiments from 2023 under various operation
scenarios. One of the underlying operation scenarios is the high-beta operation under full noninductive current drive conditions (steady-state high-beta operation) [1]. The operation will take placewith impurity seeding to reduce the heat load towards the divertors and with keeping the separatrix
density low enough to match the core operation condition. On the other hand, the concentration of the
impurities in the core causes a harmful effect to sustain the high-performance plasma by radiation
cooling and dilution. Therefore, it is important to establish a method to control the impurity transport
in the core and SOL/divertor regions. So far, a favourable spatial distribution of radiation power has
been achieved by Ne + Ar (mixed gas) seeding experiment in JT-60U, i.e. the radiated power is almost
localized in the divertor region and as a consequence, better energy confinement than Ne-only and Aronly cases has been achieved [2]. However, the interpretative and predictive simulation of such mixed
gas seeding operation has not been performed. This is partly because the number of impurity species
that the previous version of integrated divertor plasma transport code SONIC [3,4] could kineticallysolve was limited to one. The issue has been essentially resolved by restructuring SONIC code with
Multiple-Program Multiple-Data (MPMD) framework [5], which allows SONIC to calculate transport
processes of two impurity species by a kinetic impurity transport code IMPMC. The extended SONIC
code was applied to the analysis of the JT-60SA divertor plasma with two impurity species, i.e., theintrinsic C and seeded Ar gas impurities, and demonstrated the radiative divertor plasma scenario [5].
The purpose of this paper is to examine the effects of different impurity seeding species in the JT60SA divertor plasmas step-by-step in order to study the potential impurity seeding operation regime.
Aiming for this purpose, the SONIC code has been further extended to handle three or more impurity
species kinetically based on the MPMD framework mentioned above. Now the SONIC code is capableof calculating the mixed seeding impurities Ne + Ar and intrinsic C transport by IMPMC. The
impurity-impurity interaction such as the physical sputtering of C by Ne and Ar bombardment has
been also implemented. The effects of seeding mixture of Ne and Ar in JT-60SA steady-state highbeta operation scenario is demonstrated by means of the extended version of SONIC. The impurity
accumulation in the core plasma and its influence on the core plasma are also evaluated by the core
integrated code TOPICS using boundary condition of impurity contents in the core edge region
evaluated in SONIC, as the same way as in ref. [6]2. Calculation conditions
Figure 1 shows the computational grid for JT60SA. The locations of the injection ports of each
impurity species are also shown in Fig. 1. The
intrinsic C impurities are also considered in thecalculation. In this study, the parametric survey of
the Ne seeding rate (0.0 (Ar-only) – 0.28 Pa m3/s) is
performed with other parameters kept the same asref. [5] by SONIC. The Ar seeding rate is
automatically adjusted to keep the total radiation
power to be 13MW in the plasma.
The TOPICS calculation conditions are the same as
Fig. 1 Computational grid for JT-60SA.
ref. [6]. The CDBM is employed as the anomalous
The locations of the injection ports of eachheat transport model. The same transport
impurity species are also indicated.
coefficients as ref. [6] are employed for the impurity
transport calculation in the core.

C radiation
Ar radiation
Ne radiation

6

4

2

(b)

6

Prad, OD [MW]1(a)Prad, ID [MW]Prad, SOL [MW]

2

(c)

4

2

Figure 2 The stacked plot of radiation power of each impurity species as a function of Ne seeding rate:(a) the SOL region above X-point, (b) the inner divertor region, and (c) the outer divertor region.

3. Results
3.1 SONIC calculation
Figure 2 shows the impurity radiation power as a function of Ne seeding rate in (a) SOL region above
X-point, Prad, SOL, (b) inner divertor region Prad, ID and (c) outer divertor region Prad, OD. The rapid
decrease of the Ar radiation power from 1.2 MW to 0.9 MW as the Ne seeding rate increases from 0.0
to 0.02 Pa m3/s is due to the change of the balance between the friction force and the thermal force as
discussed in ref. [7]. Then, the Ar radiation power in the SOL region gradually decreases as the Ne
seeding rate increases. The Ne radiation takes place mainly in the inner divertor region and Prad, ID
increases due to the Ne seeding as seen from Fig. 3 (b). On the other hand, Ne radiation power does
not so much increase in the outer divertor region as seen from Fig. 3 (c). In summary, Prad, SOL decreases
and Prad, ID and Prad, OD increases. Such compression of impurity radiation power in the mixed-impurity
seeding case has also been observed in the JT-60U mixed impurity seeding experiment [2]. The
analysis of JT-60U experiment will be performed near future.
3.2 TOPICS calculation
The plasma/impurity transport simulation in the core has been carried out by the TOPICS code using
core edge impurity densities computed by the SONIC as a boundary parameter. The following fourcases are compared;
Case A: Ar-only (0.35 Pa m3/s), total radiation 14MW,
Case B: Ar (0.25 Pa m3/s) + Ne (0.04 Pa m3/s), total radiation 13MW,
Case B: Ar (0.20 Pa m3/s) + Ne (0.10 Pa m3/s), total radiation 13MW,
Case B: Ar (0.14 Pa m3/s) + Ne (0.20 Pa m3/s), total radiation 13MW.
It should be noted that the Case A is different from the Ar-only case regarding the toral radiation powerin Sec. 3.1. The impurity densities nAr,total and nNe,total at ρ = 0, defined as a sum of impurity densities
at each charge state, the compression of Ar cAr, defined as a ratio of the divertor Ar density tot the core
Ar density, Zeff at ρ = 0, the total impurity radiation power in the core (0 < ρ < 1), Prad,core, and Te at ρ
= 0 are summarized in Table 1. Higher cAr in Cases B, C and D than Case A has been achieved. The
lower Zeff at ρ = 0 and Prad,core than Case A has been achieved in Cases B, C, and D. The higher Te at ρ
= 0 than Case A in Cases B, C, and D has also been obtained.
TABLE 1 The summary of the impurity densities nAr,total and nNe,total at ρ = 0, defined as a sum of
impurity densities at each charge state, the compression of Ar cAr, Zeff at ρ = 0, the total impurity
radiation power in the core (0 < ρ < 1) Prad,core, and Te at ρ = 0 in each case

4. Summary
The first application of SONIC to the mixed-impurity seeding operation in JT-60SA steady-statehigh-beta plasma has been carried out to study multiple impurity seeding effect. The parametric surveyof Ne seeding rate has been performed. The Ne impurities are injected to the Ar seeded plasma with
the Ne seeding rate of 0.0 - 0.28 Pa m3/s. The rapid decrease of the Ar radiation power from 1.2 MW
to 0.9 MW as the Ne seeding rate increases from 0.0 to 0.02 Pa m3/s is due to the change of the balance
between the friction force and the thermal force as discussed in ref. [7]. Then, the Ar radiation powerin the SOL region gradually decreases as the Ne seeding rate increases. The Ne radiation takes placemainly in the inner divertor region and Prad, ID increases due to the Ne seeding. On the other hand, Ne
radiation power does not so much increase in the outer divertor region. These change of radiation
power in the SOL suggests that Ne injection into the Ar-seeded plasma has a capability of compressing
the Ar impurities into the inner divertor region. Such compression of Ar impurities into the divertor
regions results in the higher impurity radiation power in the divertor regions, and lower impurity
radiation power in the SOL region in Ar + Ne seeding case, which show the similar tendency observedin JT-60U mixed-impurity seeding experiment [2].
The core plasma has been also investigated by the TOPICS code using the core edge impurity density
computed by the SONIC code as a boundary parameter. The Ne injection into the Ar-seeded plasma
results in the higher compression ratio of the divertor Ar density to the core Ar density, i.e. the Ar
impurities have been compressed in the divertor region in Ar + Ne seeding cases. The Ar + Ne seedingcases show lower Zeff at ρ = 0 and total radiation power in the core. The higher Te than the Ar-onlycase at ρ = 0 has been obtained. These possible contributors to the better energy confinement suggestthat the mixed-impurity seeding operation is one of choices in terms of the dilution and radiation in
the core plasma.
References[1] JT-60SA Research Plan
(http://www.jt60sa.org/pdfs/JT-60SA_Res_Plan.pdf)[2] N. Asakura, et al., Nucl. Fusion 49 (2009) 115010.
[3] H. Kawashima, et al., Plasma Fusion Res. 1 (2006) 031.
[4] K. Shimizu, et al., Nucl. Fusion 49 (2009) 065028.
[5] K. Hoshino, et al., Contrib. Plasma Phys. 56 (2018) 638.
[6] N. Hayashi, et al., Nucl. Fusion 58 (2018) 066001.
[7] S. Yamoto, et al., Plasma Phys. Control. Fusion, in press.

Stability of di-vacancy and interaction with inevitable
impurities in metals
Kazuhito Ohsawa1, Takeshi Toyama2, Yuji Hatano3, Masatake Yamaguchi4
1

Research Institute for Applied Mechanics, Kyushu University
2

Institute for Materials Research, Tohoku University

3

Hydrogen Isotope Research Center, University of Toyama

4Center for Computational and e-system, Japan Atomic Energy Agency

1.

Introduction

Divertor armor tiles installed in fusion reactors will be exposed to intense plasma particle irradiation.
So, the divertor armor tiles are covered by tungsten (W) or its alloys because they have excellent
properties as plasma facing materials, e.g., low hydrogen (H) solubility, high melting point, and highmechanical toughness etc. However, a large amount of H isotopes (D and T) retention has beenreported in the W materials. In particular, tritium inventory is an important subject associated with thesafety of the fusion reactors. Vacancy-type lattice defects are important trapping sites for H isotopes.
So, the stability and binding energy of vacancy cluster with H isotopes are investigated in someprevious works [1,2]. Besides, some abnormal properties have been reported associated with vacancy
clusters in the W specimen. As an example, di-vacancy in W is energetically unstable, according tofirst-principle calculations. So, small vacancy clusters were also expected to be unstable. However,the vacancy cluster growth from mono-vacancy was observed by positron annihilation lifetime
spectroscopy which is sensitive experiment for observation of such vacancy-type lattice defects inmetals. In the present work, we propose a solution to resolve the contradiction between the simulationand experiment. Some kinds of inevitable impurities, carbon (C), nitrogen (N), and Oxygen (O) arecontained in the W materials. The vacancy clusters are remarkably stabilized by the impurities. In
particular, stabilization due to the presence O is extensively investigated.
2. Methodology
Stable structures of vacancy-type lattice defects in W and their binding energies are estimated interms of first-principle calculations on the basis of density functional theory. The simulations are
performed by using Vienna ab-initio simulation package (VASP) with PBE type potential. The cut-offenergy of plane wave is 520eV. We use large supercell composed of 432 lattice points ( 6 × 6 × 6 BCC
lattice) in the present calculations in order to reduce the effects of periodic boundary condition
imposed on the supercell. Lattice relaxations are iterated until resultant force of each atom is lowered
than 0.003 eV/Å.

The binding energies Eb of the vacancy-type lattice defects are calculated, as bellows.

Eb = E [reactants] − E [products] ,(1)

where E is the cohesive energy of the supercell. As an example, binding energy of vacancy andvacancy-impurity complex is expressed as

E b = E [Wn −1 V ] + E [Wn −1 VX ] − E [Wn − 2 V2 X ] − E [Wn ],(2)where V and V2 mean mono-vacancy and di-vacancy, respectively. Impurities are denoted by X. The
total number of lattice points for W atoms is n=432. Positive and negative binding energies correspondto attractive and repulsive interaction, respectively.
The migration energies for the impurities and mono-vacancy in the bulk W are estimated by nudged
elastic band method which is established scheme to calculate activation energy.
3. Results
We investigated two types of di-vacancy, first nearest neighbor (1nn) and second nearest neighbor(2nn) configurations, as shown in Fig. 1. Impurities in the bulk W are located close to octahedral
interstitial sites (O-sites) nest to the di-vacancy. The binding energies of di-vacancy with and without
impurities are summarized in Table I. The intrinsic di-vacancy without impurities in the bulk W is
unstable, which is good agreement with previous works [1-2]. In particular, strong repulsion isobserved at the di-vacancy in the 2nn configuration. On the contrary, the di-vacancy is remarkably
stabilized by impurities occupying an appropriate O-site. The most stable structure of di-vacancy is
the 2nn configuration and impurities are located at O-site ③, according to Table I. All impurities
contribute the stabilization but the effect of H is small, compared with other impurities (C, N, and O).
Fig. 2 shows barrier height for each impurity and mono-vacancy along the migration path in bulk
W. The migration energy for O is very small, compared with other impurities. So, O is expected toquickly diffuse in the bulk W and be trapped in a mono-vacancy at low temperature. Therefore, we
extensively investigate the effect of O for the stabilization in the present work.

Fig. 1: Schematic view of di-vacancy in 1nn and 2nn

Table I: Binding energies of di-vacancy (eV)configurations in BCC lattice. Symbols ①, ②, and ③

estimated in Eq. (2). C, N, O, and H denoteindicate O-sites next to the di-vacancy which arefavourable positions for impurity in di-vacancy.

impurities and the lowest row V corresponds tointrinsic vacancy without impurities. Symbols ①,②, and ③ indicate O-sites represented in Fig. 1.

We estimate binding energies for vacancy clusters composed of six or less than six mono-vacancies
with and without O referring to Eq. (1), as shown in Fig. 3. Small vacancy clusters without O trapping,V2 and V3, are supposed to be unstable because of the negative binding energies. However, vacancy
clusters without O composed of four or more than four mono-vacancies are stable. On the other hand,the vacancy clusters are extremely stabilized by O trapping.

4. Summary and discussion
The di-vacancy in the bulk W was expected to be unstable, according to first-principle calculations
but vacancy-cluster growth was observed in positron annihilation lifetime spectroscopy. We focus on
the inevitable impurities in W specimen and investigated the stabilization effects by the impurities.
The vacancy clusters are considerably stabilized by the impurities (C, N, O, and H). Among them, O
is expected to play the most effective role for the stabilization due to its large diffusivity.

Acknowledgement
This work was supported by JSPS KAKENHI Grant Number JP17K06993 and Computational
Simulation Center of International Fusion Energy Research Center (IFERC-CSC), Aomori, Japan in
National Institutes for Quantum and Radiological Science and Technology.

Fig. 2: Barrier height for impurities and mono- Fig. 3: Binding energies for vacancy clusters withvacancy along migration path (eV).
and without O (eV).

References[1] L. Ventelon et al., J. Nucl. Mater. 425 (2012) 16.
[2] K. Ohsawa et al., J. Nucl. Mater. 527 (2019) 151825.

Tungsten transport analysis using integrated code TASK
I. Shimomura1, N. Kasuya1,2, A. Fukuyama3, S. Mochinaga1 and M. Yagi4
IGSES, Kyushu Univ.1, RIAM, Kyushu Univ.2, Kyoto Univ.3, QST4Introduction
One of the problems for maintaining high-temperature plasmas is that impurities
generated from the wall of the reactor flow into the core plasma. Since confined energy
is lost by emission from impurities, the inflow of impurities makes plasma temperature
lower. In ITER [1] and the next-generation fusion reactor, usage of materials with highmelting point, high thermal conductivity, and high charge number, such as tungsten, isplanned [2]. Since the effect of the impurities is serious, it is necessary to control its flow
and accumulation appropriately.
In this research, code development of an integrated transport simulation is carried out
to establish control scenario of impurities in core plasmas of fusion devices. The transport
calculation routine including the transport process and atomic process of tungsten is
implemented in TASK [3] integrated transport simulation code, which enables evaluationof impurity transport processes and its response to plasma heating.

NCLASS (neoclassical transport model)
NCLASS [4] is a code for calculating neoclassical transport properties ofaxisymmetric plasmas of arbitrary aspect ratio and geometry. By giving temperature anddensity profiles of electron and ion as input parameters, the coefficients for theneoclassical particle flux of impurity ion in the following equation are calculated;
𝛤

𝐷

𝜕𝑛
𝜕𝜌

𝐷 𝑛

𝑔

𝜕𝑛
𝜕𝜌

𝑔

𝜕𝑇
.
𝜕𝜌

1

The temperatures of impurities are assumed to be equal to that of bulk ions. These
coefficients are introduced into the transport code.

Atomic process database OPEN-ADAS
The TASK code is developed to include a library for ADAS database [5] for atomic
processes such as ionization and recombination. Figure 1 shows the charge number
dependence of the difference between the ionization and recombination coefficients oftungsten obtained from the OPEN-ADAS database. Using these data, time evolution of
the impurities is solved.

Fig. 1 Charge number dependence of ionization and recombination rates.

Simulation results
Using the TASK / TR and TI transport modules, radial profiles of the density of bulk
ions, electrons and impurities in the steady state are evaluated. Here TR and TI are
modules for bulk plasmas and impurities, respectively. Figure 2 shows an example ofradial profiles of the convective velocity of bulk plasma and tungsten ions. The particles
have negative convective velocities, which corresponds to existence of inward pinch, inthis case.

Fig. 2 Radial profiles of the convective velocity of bulk plasma (left) and tungsten ions (right).

Radial profiles of the tungsten density are evaluated with fixed bulk plasma profiles
using TASK / TI module (Fig. 3). Tungsten flowing into the plasma from the boundary istransported to the core accompanied with repeated ionization and recombination. The
density accumulation of highly charged ions is enhanced by the pinch effect. Figure 4
evaluates the density peaking of the impurity by the particle pinch. As shown in Fig. 4,

the higher the inward convection velocity is, the stronger the density accumulation
becomes. This dependency is evaluated with constant convective velocities in space toshow the necessary magnitude of convection for density accumulation.

Fig. 3 Radial profile of tungsten ions.

Fig. 4 Dependence of density accumulation on the pinch velocity.

References[1] ITER Organization website, https://www.iter.org/[2] R.A. Pitts, et al., J. Nucl. Mater. 438 (2013) S48.
[3] TASK code website, http://bpsi.nucleng.kyoto-u.ac.jp/task/[4] W. Houlberg, et al., Phys. Plasmas 4 (1997) 3230.
[5] ADAS website, http://www.adas.ac.uk/

統合輸送シミュレーションスキームを用いた実験室プラズマの解析
Analysis of laboratory torus plasma using the integrated simulation scheme
N. Kasuya 1,2, I. Shimomura 2, S. Tomimatsu 2, S. Mochinaga 2, A. Fukuyama 3, M. Yagi 4, Y.
Nagashima 1,2 and A. Fujisawa 1,2
1
Research Institute for Applied Mechanics, Kyushu University, Kasuga, Fukuoka 816-8580,
Japan
2
Interdisciplinary Graduate School of Engineering Sciences, Kyushu University, Kasuga,
Fukuoka 816-8580, Japan
3
Kyoto University, Nishikyo-ku, Kyoto 615-8540, Japan
4
National Institute for Quantum and Radiological Science and Technology, Obuchi,
Rokkasho, Aomori 039-3212, Japan
E-mail: kasuya@riam.kyushu-u.ac.jp
PLAsma Turbulence Observatory (PLATO) is a torus device, which is now being underconstruction in Kyushu university [1]. The aim of the research project is to deepen
understanding of plasma turbulence phenomena by measuring spatio-temporal structures ofplasma fluctuations in this device with various high resolution diagnostics. For prediction andinterpretation of plasma phenomena in PLATO, we are carrying out simulation researches byusing various modules in TASK code [2]. TASK code is an integrated code to calculate
self-consistent states of burning plasmas by combining several modules for physical processes
in plasma core and peripheral regions. Communication between the modules is possible withinterfaces for the common data set. Figure 1 shows the schematic of module configuration in
TASK used in this research. EQ module is to calculate plasma equilibrium with a fixedboundary condition. TR module is to calculate profiles of bulk plasmas considering
neoclassical and anomalous transport processes. TI module is developing in this research to
calculate profiles of impurity ions considering atomic processes between various charge states[3]. WR module calculate RF wave propagation, so can evaluate power deposition by ECH in
tokamak. In addition to the prepared function, extension is being made to simulate
experimental diagnostics for direct comparison between simulations and experiments.
As the application of prepared functions of TASK, evaluation of equilibria and transport
analyses in PLATO are carried out. Figure 2 shows an example of the equilibrium withconsidering the external coil configuration in PLATO. This is the case with toroidal magnetic
field Bt = 0.2 T, plasma current Ip = 40 kA, central beta 0 = 0.6 %. By changing the plasma
current and external coil currents, the corresponding safety factor q profiles and plasma
shapes are calculated, such as low central q case and low triangularity case. one-dimensionaltransport analysis using TR module shows central ion temperature 50 eV and electron
temperature 200eV in the case with central electron density 11019 m-3 sustained by particle
supply in the edge region. The information is also used for inputs to other codes. Global
simulations of MHD instability and micro instability in real magnetic field geometry are
effective to capture 3-D feature of plasma perturbations. MHD simulations are carried out by
MIPS code where equilibrium profiles are introduced from TASK/EQ module. Nonlinear
saturation dynamics are represented with coexistence of ballooning and kink instabilities [4].
For the new application of the integrated code, experimental diagnostics are simulated
numerically. We have been carrying out simulations of turbulence diagnostics in experimental

devices [5]. Perturbation field that is the target for the numerical diagnostics is calculatedusing reduced MHD model in the equilibrium from TASK/EQ. Three heavy ion beam probes(HIBPs) [6] are planned to observe spatio-temporal patterns of the potential and density
fluctuations inside the PLATO plasma, and their measurements are simulated (Fig. 3). For
HIBP simulation, not only the static configuration of the plasma device including theequilibrium but also dynamics of plasma turbulence are taken into account.
Three-dimensional time-series data from a four-field reduced MHD model [7] are used for thenumerical diagnostic. Nonlinear couplings between MHD modes as tearing instability anddrift wave modes, which have different spatial scales, can be calculated. The developed codeis flexible to apply to other three-dimensional turbulence simulations and device
configurations. Microwave reflectometry is also simulated by using TASK/WR module.

Fig. 1: Schematic of module configuration in TASK integrated simulation code used in thisresearch.

Fig. 2: Equilibrium in PLATO calculated by TASK/EQU module. The positions of the verticalfield coils are also indicated.

Fig. 3: Trajectory of injected heavy ions in HIBP measurement.
Acknowledgements
This work was partly supported by JSPS KAKENHI Grant Number (JP16K06938,
17H06089), by the collaboration program of NIFS (NIFS19KNST144, NIFS18KNXN373)
and RIAM of Kyushu University. JFRS-1 supercomputer in IFERC-CSC and Plasma
Simulator in NIFS were used for the nonlinear calculations.
References[1] A. Fujisawa, AIP Conf. Proc. 1993 (2018) 020011.
[2] http://bpsi.nucleng.kyoto-u.ac.jp/task/[3] I. Shimomura, et al., this meating, 2-6.
[4] S. Tomimatsu, et al., this meating, 2-2.
[5] N. Kasuya, et al., Plasma Sci. Tech. 13 (2011) 326.
[6] T. Ido, et al., Rev. Sci. Instrum. 77 (2006) 10F523.
[7] M. Yagi, et al., Plasma Fusion Res. 2 (2007) 025.

17th Burning Plasma Simulation Initiative (BPSI) Meeting

Development of transport prediction scheme basedon first-principle simulation and machine learning
Masanori Nunami
National Institute for Fusion Science, Toki 509-5292, Japan
Based on the ﬁrst-principle gyrokinetic simulations, the machine-learning techniques, andthe reduced transport model, a new numerical scheme to predict the plasma turbulent transportlevels and the plasma proﬁles is developed. In the developed scheme, we apply the machinelearning method to ﬁnd the relevant input parameters of the nonlinear gyrokinetic simulations
which should be performed via the reduced transport model, and to optimized the reduced
transport model. Using the optimized reduced model obtained from the developed scheme, we
can predict the turbulent transport ﬂuxes by ﬁrst-principle simulations as few times as possible.

1 Introduction
To predict the turbulent transport, the ﬁrst-principle simulation based on the gyrokinetics isreliable way. Indeed, it becomes possible to validate the gyrokinetic simulation results against theexperimental observations for the plasma temperature and density proﬁles with the experimental
errors taken into account[1]. In the validation analyses, by performing gyrokinetic simulations
with changing the input parameters about the proﬁles of the plasma temperature and density, we
should ﬁnd the transport ﬂuxes that agree with the ﬂuxes from the experimental observation.
In such validation scheme, which is called “ﬂux-matching technique[2],” we can predict thetransport levels or the plasma proﬁles, quantitatively. However, to obtain quantitatively matchedtransport ﬂux with the experimental results, numerous gyrokinetic simulation runs should bedemanded. In particular, for the case of multi-species plasmas, we must perform a huge numberof the simulations to ﬁnd the matched ﬂuxes in the huge multi-dimensional parameter spaces. On
the other hand, the reduced transport model, which is based on the nonlinear simulation resultsor the data analyses via the machine-learning techniques, enables us to obtain the turbulent
transport ﬂuxes without additional heavy gyrokinetic simulations. However, since such reduced
models contain ﬁnite errors to evaluate the transport ﬂuxes in principle, it is hard to predict
the turbulent transport ﬂuxes precisely.

author’s e-mail: nunami.masanori@nifs.ac.jp

17th Burning Plasma Simulation Initiative (BPSI) Meeting

2 Developed scheme
The developed new scheme to predict the turbulent transport is proceeded as follows. In orderto reduce the number of the gyrokinetic runs, we have to ﬁnd the relevant input parametersof the simulations to be able to realize the transport ﬂuxes which is close to the experimental
observations. Therefore, ﬁrst, for ﬁnding the input parameters when the nonlinear gyrokinetic
simulation should be performed, we employ the reduced ion heat transport model[3] including
the turbulent contribution L and the zonal-ﬂow contribution τZF as represented byχmodelA1 Lα0i=,(A2 + τZF /L1/2 )χGBicrwhere, L ≡ a(ρ)[R/LT i − β0 R/Lcr
T i ] with the critical temperature gradient R/LT i (ρ) obtainedby the linear stability analyses, τZF is the normalized zonal-ﬂow decay time estimated by thelinear gyrokinetic calculation, a(ρ) is determined from the linear stability analyses, and A1 , A2 ,α0 , and β0 are constant numbers in the model. Using one of the optimization methods in themachine-learning numerical library[4], the ﬂux-matching is performed in the model to guess theinitial temperature gradient which realizes the transport ﬂux that agree with the observationsin the target plasma. Next, at the guessed temperature gradient as an input parameter, one runof the nonlinear gyrokinetic simulation is performed at each radial position. At last, using theresults of the nonlinear runs, the machine-learning is performed again to optimize the reduced
model by optimizing the parameters, α0 and β0 to α and β, respectively. Then we can obtain
the optimized reduced transport model suitable for the target plasma by only one gyrokinetic
simulation at each radial position.

3 Prediction by the developed scheme
Figure 1 shows the comparison of the temperature gradient dependences of the ion heat diﬀusivities obtained by the optimized reduced model, the initial reduced model, and the nonlinear
gyrokinetic simulations by using GKV code[5] in the case of the ion temperature gradient driven
turbulent transport in the LHD[6]. While the initial reduced model does not agree with thenonlinear simulation results, the results of the optimized reduced model with α/α0 = 0.89 and
β/β0 = 0.92, which are determined by the machine-learning, are quite agree with the nonlinear
runs around marginal temperature gradients. If the ﬂux-matching is performed using the optimized model, we can obtain the matched transport ﬂux or the matched heat diﬀusivities. In
the results of χi by the ﬂux-matching technique, the turbulent ion heat diﬀusivities obtained by
optimized model can reproduce the original nonlinear simulation results. On the other hand,
the initial model cannot reproduce the nonlinear results because the ﬁnite errors that essentiallyexists in the initial reduced model. At least in this case, we can obtain the turbulent transport
ﬂuxes which quantitatively agree with the experimental results by using the developed schemewith just once nonlinear simulation at each radial position.

17th Burning Plasma Simulation Initiative (BPSI) Meeting

Fig. 1 Comparison of the temperature gradient dependences of the ion heat diﬀusivities obtained fromthe optimized model (orange dashed curve), the nonlinear ﬁrst-principle simulations (green circles),and the initial reduced model (blue dotted curve) at ρ = 0.65 in the high-Ti LHD plasma. The
green colored region represents the experimentally allowable region of the temperature gradient.

4 Summary
The scheme to predict the turbulent transport ﬂuxes is newly developed. We employ themachine-learning to ﬁnd the relevant input parameters of the nonlinear gyrokinetic simulations
which should be done. By performing just once ﬁrst-principle simulation with the guessed input
parameters for each radial position, the optimized transport model suitable for the target plasma
can be obtained. Since the developed scheme enables us to drastically reduce the numbers of
the nonlinear simulations, the optimized model can be applied to an integrated transport code.

Acknowledgment
This work is supported in part by the RIAM collaborative Research Program (2019FP-14),
JSPS KAKENHI (18H01202, 19H01879 and 19K03801), NIFS Collaborative Research Program(KNST145, KNTT045 and KNTT050), and by FLAGSHIP2020, MEXT within the priority
study 6. The results were obtained by using JFRS-1 at IFERC CSC (KML3D and GDKTHEL),and ”Plasma Simulator” at NIFS.
[1] M. Nunami, et al., Phys. Plasmas 25, 082504 (2018).
[2] T. Görler, et al., Phys. Plasmas 21, 122307 (2014).
[3] S. Toda, et al., J. Phys. Conf. Ser. 561, 012020 (2014).
[4] M. Abadi, et al., arXiv:1603.04467 (2016).
[5] T.-H. Watanabe and H. Sugama, Nucl. Fusion 46, 24 (2006).
[6] Y. Takeiri, et al., Nucl. Fusion 57, 102023 (2017).

Neural-network-based multi-channel turbulent transportmodeling
E. Narita1 , M. Honda1 , M. Nakata2,3 , M. Yoshida1 and N. Hayashi111National Institutes for Quantum and Radiological Science and Technology, Naka,
Ibaraki 311-0193, Japan
2
National Institute for Fusion Science, Toki, Gifu 509-5292, Japan
3
The Graduate University for Advanced Studies, Toki, Gifu 509-5292, Japan

Introduction

Since density and temperature profiles severely aﬀect fusion power, understanding and predicting profile formation mechanisms are helpful to obtain high-performance fusion plasmas. A
variety of profile shapes has been found in experiments, and we have tried to predict them withseveral transport models. While theory-based turbulent transport models, such as TGLF [1],show reasonable agreement with experiments, the existing models are hard to use to investigate
the relationship between profile formation and transport processes because they do not tell us
the breakdown of the calculated particle and heat fluxes. Furthermore, introducing the detailed
descriptions of turbulence physics to the transport models tends to increase the computationalcost. To resolve the issues, we have developed the turbulent transport model DeKANIS.

2

Prediction of particle and heat fluxes consistent with thequasilinear transport matrix

DeKANIS has originally been developed for particle transport [2], and the capability has
been extended to cover both particle and heat transport, predicting turbulent particle andheat fluxes for electrons in a quasilinear limit as Γ̄e = D̄(R/Lne + CT R/LTe + CP ) and Q̄e =χ̄e (R/LTe + CN R/Lne + CHP ), respectively. Here, non-dimensional D̄, χ̄, R/Ln and R/LT are
the particle and heat diﬀusivities in proportion to fluctuation amplitudes and the density andtemperature gradients, respectively. The fluxes are composed of the diagonal (diﬀusion) andoﬀ-diagonal (pinch) terms, which are quantitatively determined by calculating the oﬀ-diagonalterm coeﬃcients (CT , CP , CN and CHP ) and the diﬀusivities (D̄ and χ̄). DeKANIS estimatesthe oﬀ-diagonal-term coeﬃcients and the diﬀusivities with a neural network (NN) model, which
is one of the machine-learning techniques, for the speed-up.
The training dataset for the NN construction has been developed with 23 JT-60U ELMy
H-mode plasmas. The plasma current and the toroidal magnetic field are IP = 1 MA andBϕ = 2.0–2.2 T, respectively, and the plasmas are heated by neutral beams (NBs) whose power
is PNB = 8.9–14.7 MW, without gas puﬃng. Plasma parameters in 0.3 ≤ ρ ≤ 0.6 are used
for development of the training dataset, where ρ is the normalized minor radius. The oﬀdiagonal-term coeﬃcients in the training dataset have been given by linear calculations of thelocal flux-tube gyrokinetic code GKW [3]. As for the particle diﬀusivity, two diﬀerent saturation
rules of fluctuation amplitudes are applied as shown in section 3: a semi-empirical saturation ruleand a mixing-length-like one. Since the remaining heat diﬀusivity is intrinsically determined tosatisfy the Onsager symmetry between the oﬀ-diagonal terms of the particle and heat fluxes, we
require the NN model that predicts the oﬀ-diagonal-term coeﬃcients and the particle diﬀusivity.

3

Two diﬀerent saturation rules of fluctuation amplitudes

The first saturation rule uses
experimental turbulent particle
fluxes Γ̄e,exp , which are estimated on the assumption thatthe external particle source flux
from the NB injection is balanced with the neoclassical particle flux and the turbulent
one. The particle diﬀusivity has
Figure 1. Radial profiles of decomposed (a) particle and (b)been calculated to match Γ̄e,exp
heat fluxes predicted by DeKANIS-1 for a plasma outside
in conjunction with CT,P calcuthe NN training dataset.
lated by GKW and R/Lne ,Te estimated by measurements. The resultant particle diﬀusivities are used as the training data withthe oﬀ-diagonal-term coeﬃcients. The NN model has been constructed to reproduce the training dataset, and the heat diﬀusivity is determined not to break the Onsager symmetry. We
note that the saturation level of the fluctuation amplitude is estimated based on Γ̄e,exp , without
the use of the experimental heat flux Q̄e,exp . The semi-empirical transport model is referred toas DeKANIS-1. Figure 1 demonstrates that DeKANIS-1 gives radial profiles of diagonal andoﬀ-diagonal terms individually, and the total particle and heat fluxes are quite close to theexperimental ones as shown in figure 2.
The second saturation rule isbased on the mixing-length theory.
In this paper, the following scaling
formula [4] is used: D̄ = 1.86 ×()1.5 −2.71
10−4 γ̄/k̄θ2
LZF , where γ̄ andk̄θ are the normalized local maximum linear growth rate and thepoloidal wavenumber corresponding to γ̄, and LZF is the residFigure 2. Radial profiles of (a) particle and (b) heat fluxes
predicted by DeKANIS-1 and -2 for the same test plasma
ual zonal flow level. Here, LZF
as in figure 1.
is given as the product of linear
zonal flow response function KRH
2
and the zonal flow( potential, which) is assumed to be proportional to a square root of γ̄/k̄θ [5],where KRH = 1/ 1 + 1.6q 2 /ϵ1/2 [6], and q and ϵ are the safety factor and the inverse aspectratio, respectively. To apply the scaling formula, the other NN model has been constructed toreproduce γ̄ and k̄θ that are given by the linear calculations. This theory-based model is referred
to as DeKANIS-2. As shown in figure 2, DeKANIS-1 and -2 reproduce the particle and heat

fluxes to a similar degree.

4

Checks of the model accuracy

Both transport models give fluxes in ∼ 10−3 second with a single CPU, which run 103 ∼ 104
times faster than the conventionally-used transport models, and they have been validated, using
14 JT-60U ELMy H-mode plasmas as test data. Here, the test data consists of 13 low-IP plasmas
with IP = 1 MA and Bϕ = 1.9–2.2 T, which are close to those of the training dataset and one
high-IP plasma with IP = 1.15 MA and Bϕ = 2.6 T.
DeKANIS-1 and-2 have beenused to predict
particle and heat
fluxes of the
training and test
data.
Here,
TGLF has also
been performed
for comparison.
For each plasma,mean standard
deviations S forthe particle andheat fluxes havebeen estimated.
In figure 3, one
Figure 3. Histograms of the mean standard deviations S of the (a)–(c)can see thatparticle and (d)–(f) heat fluxes predicted by (a), (d) DeKANIS-1, (b),(e) DeKANIS-2 and (c), (f) TGLF for the 23 plasmas applied to trainSΓ̄e ,Q̄e of theing the NNs (training, shown by red bars) and the 14 plasmas outside
two models disthe training
dataset (test, shown in blue and cyan). S for ξ is given as√(
tributes in a
) (∑
)
∑
2
2 (ρ) , where ξ
(ξ
(ρ)
−
ξ
(ρ))
/
ξ
S
=
sim
exp
sim(exp) is the simuρρ exp
similar range
between the train- lated (experimental) value of ξ.
ing and test
data, and both DeKANIS versions have the higher accuracy than TGLF.

References[1] G.M. Staebler et al., Phys. Plasmas 12 102508 (2005).
[2] E. Narita et al., Nucl. Fusion 59 106018 (2019).
[3] A.G. Peeters et al., Comput. Phys. Commun. 180 2650 (2009).
[4] E. Narita et al., Plasma Phys. Controlled Fusion 60 025027 (2018).
[5] M. Nunami et al., Phys. Plasmas 20 092307 (2013).
[6] H. Sugama et al., Phys. Plasmas 13 012501 (2006).

Data assimilation system based on integrated transport
simulation of LHD plasma
Y. Morishita1, S. Murakami1, M. Yokoyama2,4, G. Ueno3,4
1Department of Nuclear Engineering, Kyoto University, Nishikyo, Kyoto, 615-8540, Japan
2 National Institute for Fusion Science, National Institutes of Natural Science, Toki 509-5292, Japan
3The Institute of Statistical Mathematics, Research Organization of Information and Systems, Tokyo 190-8562,

Japan
4The Graduate University for Advanced Studies, SOKENDAI, Toki 509-5292 and Tokyo 190-8562, Japan

1. Introduction
Integrated simulation for fusion plasmas has various uncertainties in each of theemployed simulation models, specifically turbulent transport model. Because of this, thesimulation results also have uncertainties and we can not evaluate the employed simulation
models. To solve this problem, data assimilation techniques are applied to the integrated
transport simulation by TASK3D [1] for a plasma in LHD.
The purpose of this study is to develop a numerical system that can predict the
behavior of fusion plasmas with high accuracy employing the data assimilation. Dataassimilation can also improve simulation models involving turbulence models. As data
assimilation methods, We use the ensemble Kalman filter (EnKF) [2] to predict the behaviorof fusion plasmas with high accuracy and the ensemble Kalman smoother (EnKS) [2] toestimate the models which can reproduce experimental data. The EnKS optimizes the statevariables using both past and future data, while the EnKF using only past data.

2. Data assimilation system
In TASK3D, one dimensional (1D) diffusive transport equation is solved. As aturbulent transport model for the LHD plasma, the gyro-Bohm model for electron and the
gyro-Bohm gradT model for ion are employed based on previous TASK3D simulations as Eq.
1 and 2.
T ρ
χeTB = Ce e i
(1)
eB a

Ti ρi ∇Tia
(2)eB a ( Ti )where B is the magnetic field strength, ρi is the ion Larmor radius, a is the plasma minor
radius, and Ce, Ci are constant factors that are chosen to fit the experimental result, Ce = 23.0
and Ci = 9.07.
The state vector is composed of the electron and ion temperature, density, numericalcoefficients of turbulent transport models, and heat deposition. The NBI heat deposition isevaluated by GNET-TD [3]. Every state variable is defined on 60 computational grid points(radial direction, 1-D). The time series data of experimentally measured temperature anddensity profiles are assimilated into TASK3D.
χiTB = Ci

Te

Ti

Fig. 1 : Prediction and filtered estimate of temporal change of Te and Ti for ρ = 0.0 and ρ = 0.6.
Here, ρ is the normalized minor radius. The bar graphs show the estimated standard deviationsof background error and observation noise before filtering.

3. Simulation results
We apply the data assimilation system to a typical high ion temperature NBI plasma in
LHD (shot:114053). Figure 1 shows the simulation results of Te and Ti by EnKF, that havebeen performed with 2000 ensemble members for the cycle of assimilation, τDA = 80 msec.
The assimilation system can accurately reproduce both the electron and ion temperature interms of those profiles and temporal variations through the optimization of the employedturbulent transport model.
Figure 2 shows the filtered estimate of the numerical coefficient for electron turbulent
transport model, Ce . The radial profile have been optimized to reproduce the experimental
temperature data, but this estimate is based on only the past data. The EnKS corrects thefiltered estimate as Fig. 2 using future data. The smoothed estimate of Ce by the EnKS isshown in Fig. 3. It seems that the position of peak moves in the direction of time axis (form
1.2 sec to 1.0sec). Figure 4 is the simulation results of Te and Ti by TASK3D using thesmoothed estimates of Ce , Ci and the NBI heat deposition. It can be confirmed that theTASK3D simulation using this smoothed estimates can reproduce the experimental
temperature data with high accuracy. If the estimates of C for various plasmas can be
reproduced by a parametric or nonparametric model, more accurate turbulent transport model
can be obtained with relevant physical interpretation.

Ce /Ce0
1.4

1.0
1.0

0.0

0.5

ρ = r /a0.9

Fig. 2 : Filtered estimate of Ce by EnKF,
which are divided by the conventional value,
Ce0=23.0.

Ce /Ce01.6Time [sec]Time [sec]

1.6

1.4

1.0
1.0

0.0

0.5

ρ = r /a0.9

Fig. 3 : Smoothed estimate of Ce by EnKS,
which are divided by the conventional value,
Ce0=23.0.

Fig. 4 : Simulation results of temporal change of Te and Ti for ρ = 0.0 and ρ = 0.6 using
the smoothed estimate of Ce, Ci and the NBI heat deposition.

4. Summary
Data assimilation techniques have been applied to the integrated transport simulation(TASK3D) of a plasma in LHD. We use the EnKF and the EnKS as a data assimilation
methods for the estimation of state variables. The time series data of experimentally measured
temperature and density profiles have been assimilated into TASK3D. The assimilation
system has accurately reproduced both the electron and ion temperature through the
optimization of the employed turbulent transport model. These results indicate theeffectiveness and validity of the data assimilation approach for accurate prediction of the
behavior of fusion plasmas and the possibility of advanced turbulence modeling.

References[1] S. Murakami, et al., Plasma Phys. Control. Fusion 57, 119601 (2015).
[2] G. Evensen, Ocean Dynamics 53, 343-367 (2003).
[3] H. Yamaguchi, et al., Plasma Fusion Res. 8, 2403099 (2013.

TAE and Alfven-Slow Eigenmodes in Burning Tokamak Plasmas
C. Z. Cheng
Department of Advanced Energy, University of Tokyo, Japan
Talk presented at BPSI17, held at Kyushu University, December 15-16, 2019
In present tokamak plasma fusion experiments, the auxiliary heating such as
neutral beam injection, radio frequency wave heating, etc. are the main heating sources
to heat the ions and sustain the plasma steady state operations to counter the energy loss
due to anomalous energy transport. However, in burning tokamak plasma reactors, theauxiliary heating methods are insufficient to heat ions and sustain the burning plasma
operation. Instead, to maintain plasma in fusion burning condition, ions must be heatedby alpha particles, which are produced in the deuterium-tritium (D-T) fusion reaction
to maintain plasma temperature ideal for fusion reaction. Presently scientists expect thatalpha-particles to be the dominant heating source with alpha heating power P (=W/s)larger than the loss power Ploss (=Wtot/E), where W is the alpha particle energy, s isthe alpha particle slowing-down time, Wtot is the total plasma energy, and E is theenergy confinement time. Thus, alpha particles control the thermal plasma profiles, theglobal plasma stability and the energy and particle confinement. Thus, significant alphaparticle loss can quench DT burning, degrades heating and current drive efficiency, andcauses localized damage on first wall of fusion reactors. Thus, an important scientificfocus of fusion energy research has been the confinement of fast ions and alpha particles.
In mid-1980s the Toroidal Alfven Eigenmodes (TAE) [1, 2] was discovered
theoretically. A brief description of TAE frequency gap and existence of TAE is given
in Appendix A. Then, in early 1990s the excitation of TAEs by fast ions was observed
in TFTR and DIII-D tokamaks. Since then, many other types of Alfven eigenmodes
such as the Reversed Shear Alfven Eigenmodes (RSAE) [3-6] have been discovered.
The interaction of fast ions with Alfven eigenmodes have been very actively studied
and significant loss of fast ions caused by unstable Alfven eigenmodes has been
observed in almost all major tokamaks such as TFTR, JT-60U, JET, etc. Thus, inburning plasma fusion reactors, the Alfven eigenmodes are also expected to interact
with alpha particles and can possibly cause significant alpha loss from the reactor core.
Early D-T experiments in TFTR had identified unstable TAE and RSAE, but with weakgrowth rate due to small tritium fraction [7, 8]. Thus, no significant alpha loss wasreported. However, in burning plasma fusion reactors, large amount of tritium will be
used to produce significant fusion power. Thus, fast ion interaction with Alfven
eigenmodes will be a critical research area of burning plasma physics.

Tokamak experiments in high-β plasmas have also revealed new discrete
frequency modes with frequencies below the TAE frequency continuum gap and they
are called the Beta induced Alfven Eigenmodes (BAE) [9, 10]. One theory for the BAE
is the creation of the beta-induced Alfven gap, which is just below the lower continuum
that defines the TAE gap lower boundary and is obtained by employing the slow mode
approximation of the MHD model [11, 12]. However, the slow mode approximationneglects the slow mode propagation physics and is thus not appropriate. Thus, many
works were performed to study the frequency gaps due to coupling between the Alfven
wave and slow mode continuous spectra based on the MHD theory and the existenceof discrete frequency eigenmodes with frequencies below the TAE gap. Recently wehave provided a comprehensive theory of the formation of Alfven-Slow (AS) frequency
gaps and the discovery of many Alfven-Slow Eigenmodes (ASE) in Tokamaks [13].
The previously identified Beta-induced Alfven Acoustic Eigenmode (BAAE) [14] is
one type of several different discrete ASEs. The existence of ASEs is ubiquitous fornormal and reserve safety factor profiles, broad range of plasma β values, diﬀerent
plasma shapes and many diﬀerent toroidal mode numbers. The newly discovered ASEs
provide explanation to the experimentally observed BAE destabilized by fast ions withfrequencies below the TAE frequency gap. A simplified explanation of AS gaps and
ASEs is given in Appendix B.
With the theoretical discovery of the TAE and ASE, we expect these eigenmodes
will be critical for understanding the confinement of alpha particle and fast ions inburning plasma fusion reactors. TAEs were previously studied extensively, but ASEs
are just started to be studied. Because ASEs have lower frequency than TAEs, we expect
ASEs will interact with lower energy fast ions such as the slowed-down alpha particles
called Helium ash and the thermal ions more easily. In the next few years, theoreticaland experimental efforts should be addressed to study interaction of fast ions and alphaparticles with ASEs and TAEs.
Appendix A: TAE gap and TAEs
The TAEs are created due to the coupling of the shear Alfven continuous frequency
spectrum of neighboring poloidal harmonics such as m and m+1 around theq=(m+1/2)/n region, where q is the safety factor, n is the toroidal mode number and mis the poloidal harmonic number [1, 2]. The coupling forms a frequency gap, called
TAE gap, in the Alfven continuous spectrum with frequency as continuous functions of
the plasma minor radius and allows global TAEs to exist with frequencies inside theTAE gap. Figure 1 shows (a) frequency (/)2 versus plasma minor radius of an n=1

mode continuous spectrum and (b) the radius structure of poloidal harmonics of an n=1
TAE mode for a circular tokamak with q(0) = 1.048 and q(a)=2.3, where  =VA(0)/q(0)R is a normalized frequency [2]. Because TAE frequency does not intersectwith the continuous frequency spectrum, the poloidal harmonic structures are regularfunctions of radius and thus TAEs do not suffer continuum damping due to phasemixing effect. Thus, TAEs can be destabilized by the free energy of fast ion pressure
non-uniformity if the free energy can overcome the Landau damping effects of thermal
ions and electrons.

FIG 1 (a) Frequency (/)2 versus plasma minor radius of an n = 1 mode
continuous spectrum, and (b) the radius structure of poloidal harmonics of an n = 1 TAE
mode for a circular cross-section tokamak with q(0) = 1.048 and q(1)=2.3. The
frequency is normalized by  = VA(0)/q(1)R, VA is the Alfven velocity, and R is the
major radius.
Appendix B: Alfven-Slow Frequency Gaps and Alfven-Slow Eigenmodes
The Alfven-Slow (AS) frequency gaps are produced due to the coupling of theshear Alfven continuous frequency spectrum of m poloidal harmonic and the slow mode
continuous spectrum of different m poloidal harmonics [13]. For example, the couplingof the Alfven m poloidal harmonic and the slow mode m±1 harmonics causes the Alfvenm harmonic continuous spectrum to be uplifted from the uncoupled Alfven spectrum
by an amount related to the plasma pressure and the geodesic magnetic field curvature.
However, the frequency of one of the slow mode m±1 harmonics is reduced by the samepressure-curvature coupling eﬀect. but the frequency of the other slow mode harmonicis relatively unaﬀected. As a result the Alfven-Slow mode continuous spectrum isbroken up with two Alfven-Slow (AS) gaps below the TAE continuum gap are shown
schematically in Fig. 2. Coupling of the Alfven m harmonic with the slow mode m±2
harmonics also creates additional AS gaps if their frequencies cross each other. The ASgap widths increase with the plasma pressure. The creation of AS gaps allow theexistence of several new types of Alfven-Slow Eigenmodes (ASE). The radial structureof these new ASEs do not intersect with the continuous spectrum and thus do not suﬀer
continuum damping. An example of the continuous spectrum due to Alfven-slow mode
coupling showing the AS gaps and the existence of six ASEs with discrete frequencies
are shown in Fig. 3 for a tokamak equilibrium with β(0) = 5%, aspect ratio R/a = 3, q(0)= 1.0, and q(1)=2.42. Note that the discrete modes with frequencies above 200 kHz are
TAEs and RSAEs.

FIG. 2: The continuous spectrum (ω versus q) due to 3-harmonic coupling of the
Alfven m harmonic and the slow mode m±1 harmonics is shown schematically. The
dashed lines indicate the uncoupled slow mode frequency spectrum. The blue curvesindicate the uplifted Alfven m harmonic continuous spectrum. The green and upper red
curves show the coupled slow mode m±1 harmonic continuous spectrum. The lower
red curves indicate the slow mode m harmonic continuous spectrum.

FIG. 3: The continuous spectrum frequency showing AS gaps and the discrete
frequencies of ASEs for n = 5 mode versus minor radius for a tokamak equilibrium with
β(0) = 5%, aspect ratio R/a = 3, q(0) = 1.0, and q(1)=2.42. The red horizontal bars show
the frequencies of these discrete ASEs and the amplitude half width in the radial
direction with dot indicating the peaking amplitude location. The frequency is
calculated by choosing the normalized Alfven frequency fA0 = B0/2πρ(0)q(1)R =214.445 kHz.
References[1] C. Z. Cheng and M. S. Chance, Phys. Fluids 29, 3695 (1986).
[2] C. Z. Cheng, L. Chen, and M. S. Chance, Annals Phys. 161, 21 (1985).
[3] H. Kimura, Y. Kusama, M. Saigusa, G.J. Kramer, K. Tobita, M. Nemoto, T. Kondoh,
T. Nishitani, O. DA Costa, T. Ozeki, et al., Nucl. Fusion 9, 1303 (1998).
[4] B. N. Breizman, H. L. Berk, M. S. Pekker, S. D. Pinches, and S. E. Sharapov, Phys.
Plasmas 10, 3649 (2003).

[5] G. J. Kramer, N. N. Gorelenkov, R Nazikian, and C Z Cheng, Plasma Phys. Control.
Fusion 46, L23 (2004).
[6] M. Takechi, A. Fukuyama, M. Ishikawa, C. Z. Cheng, K. Shinohara, T. Ozeki, Y.
Kusama, S. Takeji, T. Fujita, T. Oikawa, et al., Phys. Plasmas 12, 082509 (2005).
[7] R. Nazikian, et al., Phys. Rev. Lett. 78, 2976 (1997).
[8] R. Nazikian, G. J. Kramer, C. Z. Cheng, N. N. Gorelenkov, H. L. Berk, and S. E.
Sharapov, Phys. Rev. Lett. 91, 125003 (2003).
[9] W. W. Heidbrink, E. Ruskov, E. M. Carolipio, J. Fang, M. A. van Zee eland, and R.
A. James, Phys. Plasmas 6, 1147 (1999).
[10] W. W. Heidbrink, Phys. Plasmas 15, 055501 (2008).
[11] A. D. Turnbull, E. J. Strait, W. W. Heidbrink, M. S. Chu, H. H. Duong, J. M. Greene,
L. L. Lao, T. S. Taylor, and S. J. Thompson, Phys. Fluids B 5, 2546 (1993).
[12] M. S. Chu, J. M. Greene, L. L. Lao, A. D. Turnbull, and M. S. Chance, Phys. Fluids
4, 3713 (1992).
[13] C. Z. Cheng, G. J. Kramer, M. Podesta, R. Nazikian, Phys. Plasmas 26, 082508(2019)[14] N. N. Gorelenkov, H. L. Berk, N. A. Crocker, E. D. Fredrickson, S. Kaye, S.
Kubota, H. Park, W. Peebles, S. A. Sabbagh, S. E. Sharapov, et al., Plasma Phys. Control.
Fusion 49, B371 (2007).

Extraction of spatial structures of intermittent eventsusing dynamical mode decomposition
1,2

M. Sasaki, 3Y. Kawachi, 2,4,5R. O. Dendy, 6H. Arakawa, 1,2N. Kasuya,
7
F. Kin, 1K. Yamasaki, 1,2S. Inagaki
1

Research Institute for Applied Mechanics, Kyushu University, Japan
Research Center for Plasma Turbulences, Kyushu University, Japan
3
Interdisciplinary Graduate School of Engineering Sciences, Kyushu University, Japan
4
Department of Physics, Warwick University, Coventry UK
5CCFE, Culham Science Centre, Oxford UK
6
Institute of Science and Engineering, Academic Assembly, Shimane University, Japan
7
National Institutes for Quantum and Radiological Science and Technology, Japan
2

Intermittent events in turbulence are ubiquitous in fluids; the solar flare [1], thetorrential heavy rain [2], ocean current, and so on. In magnetized plasmas, intermittent
events frequently occurs, such as the edge localized modes (ELMs), and blobs, whichsignificantly affect the properties of plasmas. These intermittent events are caused
nonlinearly, so that the spatio-temporal properties usually involve higher harmonics.
Therefore, the number of degrees of freedom that play roles in intermittent events
increases when the Fourier mode decomposition is applied. In order to understand thenonlinear properties of such events, it is important to extract spatio-temporal structureswith reducing the degrees of freedom as much as possible.
In this study, the dynamical mode decomposition (DMD) is applied toturbulence simulation data [3], in order to extract the spatio-temporal behavior of theintermittent or transient events with reducing the number of degrees of freedom. The
DMD method is used, which does not assume any functional form for the structure.
Only based on observed signals, spatio-temporal patterns that are statistically valid arederived from the eigenvalue problem. The turbulence simulation we use is a direct
numerical simulation in cylindrical magnetized plasmas, and the nonlinear dynamics
of the Kelvin-Helmholtz instability are self-consistently calculated. Using the
turbulence simulation, where all the three dimensional turbulence informations are
known, the intermittent or transient dynamics of the turbulence is decomposed by the
DMD method. The obtained DMD modes comprise rotationally symmetricdeformations together with spiral patterns. A new method to calculate the time
evolution of DMD mode amplitudes is proposed, based on convolution-type
correlation integrals, and then applied to the simulation outputs in a limit cycle regime.
The resulting time traces capture the essential physics far better than Fourier
techniques applied to the same data [4].
[1] N. Scafetta, et. al., Phys. Rev. Lett. 90, 248701 (2003).
[2] N. Takahashi, et. al., Journal of Disaster Research, 14, 235 (2019).
[3] M. Sasaki, et. al., Phys. Plasmas, 26, 042305 (2019).
[4] M. Sasaki, et. al., Plasma Phys. Control. Fusion, 61, 112001 (2019).

Statistical induction of a thermal transport model of fusion plasmas
based on the transport analysis database
Masayuki Yokoyama1,21National Institute for Fusion Science, National Institutes of Natural Sciences,
2
The Graduate University for Advanced Studies (SOKENDAI),
322-6 Oroshi, Toki 509-5292, Japan
yokoyama@nifs.ac.jp

A statistical induction of a thermal transport model of fusion plasmas has been conductedby regression analysis utilizing the accumulated experimental transport analysis database [1].
This approach is unique and simple in a sense that it does not require any detailed physics
modelling such as neoclassical and turbulent transport. It has come to the stage to elucidate
the influencial physics quantities to prescribe the heat diffusivity, and to compare/interprete
experimentally identified trends. This progress will foster acquiring practical transport models
for a real-time control and the guidance to the parameter dependence to be elucidated throughlarge-scale cutting-edge simulations.
The thermal transport analysis database considered in this study is the same as that used inRef. [2], as reproduced in Figure 1. This is the database (about 3000 data) for the ion thermal
diffusivity, , constructed from 31 discharges (so-called high-ion-temperature (Ti) scenario)in Large Helical Device (LHD) [3]. The database was constructed by the integrated transport
analysis suite, TASK3D-a [4] to LHD plasmas.
Fig. 1: All the data points of
employed database (LHDhigh-Ti (ion temperature)scenario) for (a) Ti, (b) Te(electron temperature), (c)ne (electron density) andthen (d) , which are
reproduced from Ref. [2].

The log(10)-linear multivariate regression analysis was performed on this database. When it
comes to a physics interpretation of the regression result, keeping a smaller number ofvariables with larger influence is more transparent. For this, the database for reff/a99<0.35 (404data, roughly inside ion internal transport barrier [5]) is employed in the following discussion.
Here, a99 is a measure of the experimentally known last closed flux surface, defined by theeffective minor radius which encloses 99% of the total electron pressure (requiring precalculated VMEC equilibrium database, and reff is the minor radius defined by the radius of
the equivalent simple torus which encloses the same volume as the flux surface of interest [6].
The AIC (Akaike’s Information Criterion) [7] is the measure for the relative quality ofstatistical models for a given dataset. Practically, the minimum (or the lowest) AIC selects an“optimal” model among many candidates. More precisely, AICc (AIC with a correction forsmall sample sizes) [8] is used. Based on the AICc analysis, the regression expression with 4
variables (selected out of initially considered 9 variables [1]),is obtained, where, R/LTi and are the ion collision frequency normalized by its bananaplateau boundary, the scale length of the ion temperature gradient with R being the vacuum
magnetic axis position and the main helicity of magnetic configuration, respectively. It should
be noted that this expression is only valid for dataset of reff/a99<0.35 of the database shown in
Fig. 1, due to the inherent nature of regression analysis.
Positive value of power of ρi*indicates gets larger as Ti becomes higher. However, powerof Ti becomes negative if ρi* and Te/Ti considered simultaneously. The itself currently does
not cover any of changes of magnetic configurations, however, its exclusion worsens the
goodness. The negative power of R/LTi would be physically interpreted as that the ion heatdiffusivity in ion internal transport barrier decreases as the ion temperature gradient becomeslarger (or its scale length becomes smaller making R/LTi larger) [5]. The increase of heatdiffusivity (then the decrease of the ion temperature gradient) for increasing Te/Ti has also
been recognized experimentally in LHD [9].
Based on these considerations, the induced regression expression statistically resolves
complicated entanglements among Ti, Te/Ti and its gradient (R/LTi), and the feature oftransport barrier formation, in a single line expression.
It would be interesting to consider this kind of regression results to compare withsimulations when results for substantial number of time slices in LHD are accumulated. Also,transport analysis datasets in other operational scenario/other experiments (worldwide) could
provide a new insight for thermal transport modelling of fusion plasmas, complementing
conventional global scalings for the energy confinement time.
References[1] M. Yokoyama. Nucl. Fusion 59 (2019) 094004. (Corrigendum: Nucl. Fusion 60 (2020)
029501.
[2] M. Yokoyama et al., Plasma Fusion Res. 9 (2014) 1302137.
[3] Y. Takeiri et al., Nucl. Fusion 57 (2017) 102023.
[4] M. Yokoyama et al., Nucl. Fusion 57 (2017) 126016.
[5] K. Nagaoka et al., Nucl. Fusion 59 (2019) 100602.
[6] C. Suzuki et al., Plasma Phys. Control. Fusion 55 (2013) 014016.
[7] H. Akaike, "Information theory and an extension of the maximum likelihood principle",
Proceedings of the 2nd International Symposium on Information Theory, Petrov, B. N., andCaski, F. (eds.), Akadimiai Kiado, Budapest (1973) 267.
[8] N. Sugiura, Communications in Statistics - Theory and Methods 7 (1978) 13.
[9] H. Takahashi et al., Nucl. Fusion 57 (2017) 086029.