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

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

および

オンライン

1 月 5 日(木)(核融合エネルギーフォーラムサブクラスターとの合同会合)
9:30 – 9:40 林（QST）
連絡事項
9:40 – 10:10 宮戸（QST）
IFERC 計算機シミュレーションセンターの現状報告
10:10 – 10:40 横山（核融合研）、林（QST）
ITER に於ける統合モデリング活動の報告
10:40 – 11:10 若月（QST）
ITPA 統合運転シナリオグループ活動報告
11:10 – 11:30 星野（慶應大）炉心プラズマ＋炉工学炉材料モデリングサブクラスター合同会合の報告
11:30 – 12:00 林（QST）
今後の予定の議論
12:00 – 13:30 昼休み
ワークショップ「不純物輸送」 Lunch time workshop “Impurity transport”
13:00 – 13:30 林 (QST)、藤田 (名大)、糟谷 (九大) online
昨年の企画セッションを受けた共同研究作業会を昼休みに開催します(関係者による集まりですので、自分は関係者だと思われる方ご参加ください)13:30 – 13:40 事務連絡 Business announcement(座長：村上)13:40 – 14:05 講演 1-1

成田 (QST)
1

Improvement in a semi-empirical turbulent transport model for its versatility
and its application to integrated simulations
14:05 – 14:30 講演 1-2
福山 (京大)Progress of the development of TASK code and interface to IMAS
14:30 – 14:55 講演 1-3
横山 (核融合研)
Introduction to IAEA coordinated research program: AI for Fusion
14:55 – 15:15 講演 1-4
森下 (京大)Control of LHD plasma using the data assimilation system ASTI
15:15 – 15:25 休憩(座長：横山)15:25 – 15:50 講演 1-5
藤田 (名大) online
Study on operation scenario with enhanced radiation power from main plasma
by impurity injection in DEMO
15:50 – 16:10 講演 1-6
持永 (九大)Impurity transport simulation in tokamak plasmas using the integrated code
TASK
16:10 – 16:30 講演 1-7
Park (SNU)
Integrated modelling of high-temperature fusion plasma regime facilitated by
fast ions in the KSTAR tokamak
16:30 – 16:45 講演 1-8
竹本 (名大) online
Validation of tungsten transport model by simulating soft x-ray measurementin JT-60U
16:45 – 16:55 休憩(座長：糟谷)16:55 – 17:20 講演 1-9
Hahm (SNU)
ExB vortex dynamics inside a magnetic island
17:20 – 17:45 講演 1-10
矢木 (QST)
Particle transport in fusion plasmas
17:45 – 18:10 講演 1-11
Choi (SNU) online
Gyrokinetic simulation of KSTAR AE-ETG system
18:10 – 18:35 講演 1-12
登田 (核融合研)
Study of turbulent transport in trapped electron and ion temperature gradientmodes in high collision regime for tokamak plasmas
18:50

散会

2

1 月 6 日(金)9:00 – 10:30 ポスター講演 Poster(online, 開室 start 8:45、閉室 end 13:00)
10:30 – 10:40 休憩(座長：村上)10:40 – 10:55 講演 1-13
麝嶋 (日大)Estimation of 2D distribution evolution of electrostatic potential fluctuations
using deep learning
10:55 – 11:20 講演 2-1
滝塚 (阪大) online
Several low-cost techniques for global plasma simulation with particle-in-cellmodel
11:20 – 11:45 講演 2-2
沼波 (核融合研)
Manifold structure appearing in plasma turbulence simulations and transportmodeling(Session Leader：村上)
11:45 – 12:05 議論 Discussion
12:05 – 12:10 事務連絡 Business announcement
12:10 – 13:20 昼休み(特定研究「計測・シミュレーション・モデリングを組み合わせた統合診断」研究集会との合同会合)(座長：糟谷)13:20 - 13:30 講演 3-1
糟谷（九大）
特定研究「計測・シミュレーション・モデリングを組み合わせた統合診断」の紹介
13:30 - 14:00 講演 3-2
糟谷（九大）
Global structure and its diagnostic simulation in magnetized plasmas
磁化プラズマにおける巨視的構造診断のためのシミュレーション
14:00 - 14:30 講演 3-3
本多（京大）
Carbon impurity modeling with the superstage model in TASK/TX
スーパーステージモデルを用いた TASK/TX における炭素不純物モデリング14:30 – 15:00 講演 3-4
佐々木（日大）Study on spatio-temporal dynamics in turbulence applying data-driven science
データ駆動科学による乱流時空間ダイナミクスの研究
15:00 – 15:30 講演 3-5
山田（九大）
Analysis of turbulent structures in magnetized cylindrical plasmas
直線磁化プラズマにおける乱流構造の解析
3

15:30 – 16:00 糟谷（九大）
discussion 議論
16:00

散会

Poster number

20th BPSI meeting (2023.1)(Ver.3)

PA-1 村上定義 (京大) S. Murakami
Effects of magnetic and electrostatic fluctuations on ECH supra-thermal electron
behavior and toroidal torque in tokamak plasma
PA-2 轟晴彦 (九大) H. Todoroki
Analysis of turbulent particle fluxes by using 5-field RMHD code in torus plasmas
PA-3 梅崎大介 (九大工) D. Umezaki
Effect of large-angle scattering transport in divertor plasmas
PA-4 矢本昌平 (量研) S. Yamoto
Improvements of impurity transport modelling by integrated divertor code SONIC and
its application to predictive simulation of future devices
PA-5 吉原稜 (九大) R. Yoshihara
Simulation of 3D perturbation measurement in torus plasmas using a heavy ion beamprobe
PA-6 古田原拓実 (日大) T. Kodahara
Analysis of particle transport in PANTA by using multi-field SVD
PA-7 田原康祐 (京大) K. Tahara
Influence of the beam-beam collisions on the energetic particle distribution in Large
Helical Device

4Improvement in a semi-empirical turbulent transport model for
its versatility and its application to integrated simulations
E. Narita1 , M. Honda2 , M. Nakata3,4 , M. Yoshida1 , N. Hayashi1 and T. Nakayama4
1

National Institutes for Quantum Science and Technology, Naka, Ibaraki 311-0193,
Japan
2
Graduate School of Engineering, Kyoto University, Nishikyo, Kyoto 615-8530, Japan
3
National Institute for Fusion Science, Toki, Gifu 509-5292, Japan
4
The Graduate University for Advanced Studies, Toki, Gifu 509-5292, Japan

1

Introduction

Neural network (NN) based turbulent transport models have been introduced to predict temperature and density proﬁles quickly. Such models successfully predict turbulent ﬂuxes 10 3−5
times faster than the conventional reduced models, such as TGLF [1,2], and they are now utilized practically. While most NN based turbulent transport models are surrogate models, which
are built to mimic the reduced models, our model named DeKANIS [3-5] is based on gyrokinetic
calculation results and experimental data.
DeKANISﬂuxes for electrons( predicts turbulent) particle and heat
) in a quasilinear limit as
(
R
R
R
R
Γ̄e = D̄ Ln + CT LT + CP and Q̄e = χ̄e CN Ln + LT + CHP , respectively. Here, noneeeedimensional D̄, χ̄, R/Ln and R/LT are the particle and heat diﬀusivities in proportion to
ﬂuctuation amplitudes and the density and temperature gradients, respectively. The ﬂuxes arecomposed of the diagonal (diﬀusive) and oﬀ-diagonal (non-diﬀusive) terms, which are quantitatively determined by calculating the oﬀ-diagonal-term coeﬃcients (CT,P,N,HP ) and the diﬀusivities (D̄ and χ̄e ). In addition to the electron ﬂuxes, DeKANIS predicts the ion heat ﬂux, whichχ̄i,eﬀχ̄e,eﬀ LRT nnei TTei . Here, χ̄eﬀ is the normalized eﬀective heat diﬀusivity. Onceis given as Q̄i = χ̄e,eﬀiχ̄e,eﬀ is obtained with Q̄e , Q̄i can be constructed through χ̄i,eﬀ /χ̄e,eﬀ . To forecast Γ̄e , Q̄e andQ̄i , the NN model and the scaling formula are utilized. The NN model was trained to evaluate
CT,P,N,HP and χ̄i,eﬀ /χ̄e,eﬀ , whose training datasets are prepared by the linear gyrokinetic calculations with GKW [6] at the normalized wavenumber k̄θ with the normalized maximum linear
growth rate γ̄max of the ion temperature gradient (ITG) mode or the ITG/trapped electronmode (TEM). The NN model also predicts the corresponding k̄θ and γ̄max , which are passed
to the scaling formula to determine D̄. The remaining coeﬃcient, χ̄e , is calculated with therelationship between D̄ and χ̄e : χ̄e = D̄(CT + 1.5)/CN .
In this paper, the NN training datasets and the scaling formula are modiﬁed for the versatilityof DeKANIS, and it is applied to integrated simulations of an ITER pre-fusion power operation
1 (PFPO-1) plasma [5].

2Expansion of training data

The training datasets of the previous NN model were based on JT-60U H-mode plasmas, whose
ion temperatures are higher than the electron ones: Te /Ti < 1, at most radial points. Also, thedatasets cover the plasma parameters in 0.30 ≤ ρ ≤ 0.65. To apply DeKANIS to predictions of
the high-performance plasmas, it needs to be applicable to the range including Te /Ti >∼ 1 sincealpha particles primarily heat electrons in burning plasmas, electron cyclotron (EC) heating isenvisaged in the PFPO-1 phase, and the collisional equipartition leads to Te /Ti ∼ 1 in highdensity plasmas. When it comes to estimating the performance including the edge pedestal
region, DeKANIS needs to give the turbulent ﬂuxes also in the outer core region, 0.65 ≤ ρ ≤ 0.85.
To this end, a subset of JET plasmas [7] has been introduced, whose Te /Ti is around or exceeds
unity, and the outermost radial point has changed from ρ = 0.65 to ρ = 0.85.

3(Modiﬁcation of the saturation model
The previous DeKANIS determined D̄ with the following scaling formula: D̄model,pre ∝)α [()0.5 ]−β. Here, the normalized quasilinear diﬀusivity γ̄/k̄θ2 characterizes turKRH γ̄/k̄θ2γ̄/k̄θ2bulent transport, and the suppression of instabilities due to zonal ﬂows is expressed with thelinear zonal ﬂow response function KRH given as KRH = 1/(1 + 1.6q 2 /ϵ0.5 ) [8] and the proxy of()0.5
the zonal ﬂow potential, γ̄/k̄θ2, where q and ϵ are the safety factor and the inverse aspect ratio, respectively. The exponents α and β were optimized by the regression analysis to reproduce
D̄exp , which is semi-empirically estimated with the electron heat ﬂux of the JT-60U plasmas,
and α = 1.32 and β = 2.93 were obtained. As displayed in ﬁgure 1(a), the previous modelshows a high mean squared logarithmic error (MSLE) of 0.322 for JT-60U regression data, andunderestimates test data given by the JET plasmas. Furthermore, the net exponent of γ̄/k̄θ2 isnegative, α − 0.5β = −0.145, which contradicts the mixing-length theory.
Since one of the questions of the previous model is that KRH was derived in the collisionless)(2 ω 2 /ν 2 0.5 ,limit, collisional decay of the zonal ﬂows is introduced with τ̄r . Here, τ̄r = νf2 /ωf2 + KRH
df
and the square of it measures the suppression of turbulence by zonal ﬂows considering its colli()
sional damping [9], where ωf = vth,i /R, νf = ωf exp −q 2 , νd = 0.67/ (τii ϵ) with the ion thermalvelocity vth,i , the major radius R and the ion-ion collision time τii . The new scaling formula uses
(

)α [

(

)0.5 ]−βτ̄r instead of KRH : D̄model,new ∝ γ̄/k̄θ2τ̄r γ̄/k̄θ2. The optimized exponents are α =20.755 and β = 1.20, and the net exponent of γ̄/k̄θ is positive, α−0.5β = 0.155. The value is small,but the new model is better than the previous one in terms of the mixing-length theory. As shownin ﬁgure 1(b), the MSLE decreases for regression data, and the tendency to underestimate D̄exp of
test data is no longer observed. The increase in the MSLE for test data is due to the data points
with a high collision frequency, namely, νii∗ >0.1. Since enough regression data were not collected in νii∗ > 0.1, themodel predicts too largeD̄ because of the too suppressed τ̄r for the data
points in question. The
Figure 1. The particle diﬀusivity estimated by the previous model
MSLE is reduced without(a) D̄model,pre and the new model (b) D̄model,new compared withD̄exp for regression and test data.
such data points.

4

An integrated simulation of an ITER PFPO-1 plasma

The modiﬁed DeKANIS is implemented in the integrated model GOTRESS+ [10] to perform
a core-pedestal coupling integrated simulation of an ITER PFPO-1 plasma with the plasma
current of 5.0MA and the toroidal magnetic ﬁeld of 1.8T. To solve heat transport in the coreregion, ρ ≤ 0.85, DeKANIS gives the turbulent ﬂuxes and Matrix Inversion method gives theneoclassical ones. The pedestal height is determined with the in-house EPED1 model, and theequilibrium and current proﬁle alignment code ACCOME is also contained. The density proﬁle isprescribed. When the plasma is heated by EC heating with 20MW, the electron and ion temperatures are obtained as shown in ﬁgure 2(a) after eight iterations. As the turbulent ﬂuxes near themagnetic axis are computed by extrapolation, the temperatures in ρ ≥ 0.25 are compared withthose predicted by the previous study using TGLF [11], which shows Te ∼ 10keV and Ti 5.5keV
at ρ = 0.3 and Te,i ∼ 3.5keV at ρ = 0.8. Our simulation predicts similar temperature proﬁles. In forecasting the turbulent ﬂuxes, DeKANIS calculates the oﬀ-diagonal term coeﬃcients,
which can be veriﬁed with GKW as shown in ﬁgure 2(b). The ﬁrst oﬀ-diagonal term, CN R/Lne ,given by DeKANIS is close to that of GKW, while the second one, CHP , is overestimated.
Since CHP estimated by
DeKANIS and GKW areboth smaller than theother two terms, it can besaid that turbulent transport is mostly caused by
the diagonal term and theoﬀ-diagonal one driven
by the density gradientand that the two terms
Figure 2. (a) The temperature and (b) the diagonal and oﬀproduce turbulent ﬂuxes
diagonal term proﬁles predicted by the integrated simulation with
with the similar level.
DeKANIS implemented in GOTRESS+ for the ITER PFPO-1
plasma. In the subﬁgure (b), the results are veriﬁed with GKW.

References[1] G.M. Staebler et al., Phys. Plasmas 12 102508 (2005).
[2] G.M. Staebler et al., Nucl. Fusion 61 116007 (2021).
[3] E. Narita et al., Nucl. Fusion 59 106018 (2019).
[4] E. Narita et al., Nucl. Fusion 61 116041 (2021).
[5] E. Narita et al., Contrib. Plasma Phys. e202200152 (2023).
[6] A.G. Peeters et al., Comput. Phys. Commun. 180 2650 (2009).
[7] The ITER 1D Modelling Working Group, Nucl. Fusion 40 1955 (2000).
[8] M.N. Rosenbluth and F.L. Hinton, Phys. Rev. Lett. 80 724 (1998).
[9] F.L. Hinton and M.H. Rosenbluth, Plasma Phys. Control. Fusion 41 A653 (1999).
[10] M. Honda et al., Nucl. Fusion 61 116029 (2021).
[11] A. Loarte et al., Nucl. Fusion 61 076012 (2021).

Introduction to IAEA coordinated research program (CRP): AI for Fusion (AI4F)
M. Yokoyama
National Institute for Fusion Science, Rokkasho Research Center
The IAEA organized the Technical Meeting on Artificial Intelligence (AI) for Nuclear Technology
and Applications (#AI4Atoms) on 25–29 October 2021, to programmatically discuss to form theroadmap for IAEA role & support in use of AI in nuclear science, technology, and applications. It
consisted of 10 working groups (WGs) along with cross-cutting groups (enabling infrastructure, andadvanced modelling and simulation methodologies). There is a “Nuclear Fusion” WG (NFWG) in 10
WGs.
Prof. Hideo Nagatomo (Osaka University) and I (Masayuki Yokoyama) joined this NFWG fromJapan, to launch continuously discuss to form the Coordinated Research Project (CRP). It resulted in
the launch of CRP on AI for Accelerating Fusion R&D (2022–2027) with the objective of accelerating
fusion R&D with Machine Learning (ML)/AI, through the creation of a platform and cross-community
network
forinnovation
andpartnership.
(cf.,
https://nucleus.iaea.org/sites/ai4atoms/ai4fusion/SitePages/What-is-AI-for-Fusion-.aspx?web=1)
This CRP is composed by four Work Packages (WP), as follows.
WP 1: Real-time MFE System Behaviour Prediction, Identification & Optimization Using ML/AI
Methods (MFE: Magnetic Fusion Energy)
WP 2: IFE Physics Understanding through Simulation, Theory and Experiment Using ML/AI Methods(IFE: Inertial Fusion Energy)
WP 3: Feasibility of MFE and IFE Image Database
WP 4: Community Engagement & Workforce Development
The objectives, expected activities/output/outcomes for each WP are described on the abovementioned webpage.
This information was distributed via the mailing list of The Japan Society of Plasma Science andNuclear Fusion Research, for calling the proposal from institutions/universities in Japan. Then, Osaka
University and NIFS individually applied, and they were approved by the IAEA.
Below, only contents of the application from NIFS are introduced (referring the application).

Title of the research agreement
Radiation-collapse discharge database from Large Helical Device (LHD)
Summary of the proposed research
Radiation-collapse discharge database will be provided from Large Helical Device (LHD), as a part ofthis internationally coordinated research program. This should accelerate fusion R&D via ML/AI
applications regardless the differences on confinement concepts in MFE. In particular, provided
database should be a common ground for developing anomaly detection methods via ML/AI.
Description of research objectives and anticipated outcomesTo accelerate fusion R&D by establishing a multi-machine database of experimental MFE data(adhering to FAIR/Open Science principles) for ML/AI-driven applications, and through increased
access to knowledge and information of ML/AI methods for MFE. Radiation-collapse discharge
database could be provided from Large Helical Device (LHD). This should accelerate fusion R&D via
ML/AI applications regardless the differences on confinement concepts in MFE.
Scientific scope of the project
Concrete contribution to WP1 by providing radiation-collapse discharge database of LHD, as a part ofmulti-machine database of experimental MFE data for ML/AI applications.
Expected Outputs* Accelerate fusion R&D via ML/AI application.
* Ensured data reproducibility and accessibility for across domains for the international community.
* Increase access to knowledge and information about ML/AI for MFE R&D.
The progress and on-going activities within this CRP framework will be reported as they become
available.

Impurity transport simulation in tokamak plasmas
using the integrated code TASKS. Mochinaga1), N. Kasuya1,2), A. Fukuyama3) and M. Yagi4)
1)

1.

IGES, Kyushu University, 2) RIAM, Kyushu University, 3 )Kyoto University, 4)QST

Introduction

Negative magnetic shear mode and high βp mode are improved confinement modes, which can achieve highplasma performance with the formation of an internal transport barrier (ITB) in the core plasma. However, the ITB
region can increase accumulation of impurities in the core plasma, which may lead to plasma collapse [1,2]. Tungsten
produced from metal walls and divertor plates induce particularly large radiation losses, so the accumulation in thecore plasma should be avoided. One of the factors that accelerates the tungsten accumulation is inward neoclassicalpinch component driven by bulk ion density gradient [3]. Some studies have shown that electron cyclotron heating
suppresses inward impurity transport by increasing the bulk temperature gradient and decreasing the density gradient[4]. In this study, simulations are performed by focusing on neoclassical transport, which may accelerate impurity
accumulation in the ITB region, and the mechanism of tungsten profile formation is investigated.
.

2.

Transport simulation models

The TASK code is one of the integrated codes developed in Japan, which has a modular structure combining
many codes such as equilibrium, transport, and wave analysis. The multiple modules exchange data through the datainterface BPSD [5], which enables self-consistent simulations of time evolutions of plasma variables.
In the transport calculations of plasma main component using the TASK/TR code, the following diffusive onedimensional equations for particle and heat transport are solved [6].
1 ∂1 𝜕(𝑉 ′ 𝛤𝑆 )
′)
(𝑛
(1)
𝑉
=
−
+ 𝑆𝑠 ,
𝑉 ′ ∂𝑡 𝑠
𝑉 ′ 𝜕𝜌
1 ∂ 3
1 ∂(𝑉 ′ 𝑄𝑆 )
′ 5/3
(2)
(
𝑛
𝑇
𝑉
)
=
−
+ 𝑃𝑠 ,𝑠𝑠
𝑉 ′ ∂ρ
𝑉 ′ 5/3 ∂𝑡 2
Here, ns is the plasma density, Ts is the plasma temperature, Ss is the particle source, Ps is the heat source, ρ is thenormalized radius, V’ is the radial derivative of plasma volume, and s is the index representing plasma species as ion
and electron. The particle flux Γs and heat flux Qs are given as follows;𝜕𝑛𝑠(3)𝛤𝑠 = 〈|𝛻𝜌|〉𝑉𝑠 𝑛𝑠 − 〈|𝛻𝜌|2 〉𝐷𝑠
,
𝜕𝜌
and
3
𝜕𝑇𝑠 3
(4)
𝑄𝑠 = ⟨|𝛻𝜌|⟩ 𝑛𝑠 𝑇𝑠 𝑉𝐸𝑠 − ⟨|𝛻𝜌|2 ⟩𝑛𝑠 𝜒𝑠+ 𝛤𝑇 .
2
𝜕𝜌 2 𝑆 𝑠
2
Here, ˂|∇ρ|˃, ˂|∇ρ| ˃ are the metric values, Vs is the particle convective velocity, Ds is the particle diffusion
coefficient, VEs is the heat convective velocity, and χs is the heat diffusion coefficient. The diffusion coefficients are
given by linear summation of neoclassical and turbulent diffusion coefficients, and the convective velocities consistof the neoclassical component only. For the neoclassical transport model, NCLASS [7] is applied. For the turbulent
transport model, CDBM (Current Diffusive Ballooning Mode) or mBgB (mixed Bohm and gyro Bohm) models areapplied. The diffusion coefficients by the CDBM model are as the followings [8,9];𝑐 2 𝑣𝐴𝜒𝐶𝐷𝐵𝑀 = 𝜒𝑒 = 𝜒𝑖 = 𝐹(𝑠, α)α1.5 2,(5)
ω𝑝𝑒 𝑞𝑅
1(6)
𝐷𝑒 = 𝐷𝑖 = χ𝐶𝐷𝐵𝑀 ,2
where, 𝑐 is the speed of light, 𝜔𝑝𝑒 is the electron plasma frequency, 𝑣𝐴 is the Alfvén velocity, 𝑞 is the safety factor, 𝑠 ≡(𝑟⁄𝑞 )(𝑑𝑞 ⁄𝑑𝑟 ) is the magnetic shear, 𝛼 ≡ −(2𝜇0 𝑞2 𝑅⁄𝐵2 )(𝑑𝑝⁄𝑑𝑟 ) is the normalized pressure gradient. The shaping factor
𝐹(𝑠, 𝛼) is defined as

𝐹(𝑠, α) =1√2(1 − 2𝑠 ′ )(1 − 2𝑠 ′ + 3𝑠 ′ 2 )(𝑠 ′ = 𝑠 − α < 0),(7)

1 + 9√2𝑠 ′5/2(𝑠 ′ = 𝑠 − α ≥ 0).
√2(1 − 2𝑠 ′ + 3𝑠 ′ 2 + 2𝑠 ′ 3 )𝐹(𝑠, 𝛼) decreases rapidly when s-α becomes smaller than 0.75. The stabilizing effect of 𝐸 × 𝐵 flow shear is notconsidered in this study. On the other hand, the mBgB model can be expressed as follows [10,11].
(8)
𝜒𝑒 = 8.0 × 10-5 𝜒𝐵 + 3.5 × 10-2 𝜒𝑔𝐵 ,
-4
-2
(9)
𝜒𝑖 = 1.6 × 10 𝜒𝐵 + 1.75 × 10 𝜒𝑔𝐵 ,𝜒𝑒 𝜒𝑖
𝐷𝑒,𝑖 = (0.3 + 0.7𝜌).
(10)𝜒𝑒 + 𝜒𝑖
The Bohm term χB and gyro-Bohm term χgB are given by𝑎(𝑑𝑝𝑒 ⁄𝑑𝑟 ) 𝑇𝑒 (𝜌 = 0.8) − 𝑇𝑒 (𝜌 = 1.0)
1.47ω𝐸×𝐵(11)𝜒𝐵 = 𝜌𝑠 𝑐𝑠 𝑞 2
Θ (−0.14 + 𝑠 − ||),𝑝𝑒𝑇𝑒 (𝜌 = 1.0)
γ𝐼𝑇𝐺
(𝑑𝑇𝑒 ⁄𝑑𝑟 )(12)
𝜒𝑔𝐵 = 𝜌𝑠2 𝑐𝑠,
𝑇𝑒where, ρs is the gyro radius, cs is the sound speed, pe is the electron pressure, Θ is the Heaviside function, ωE×B is the
E×B flow shear rate and γITG is the ITG growth rate. In the impurity transport calculations using the TASK/TI code,
the following particle transport equations including atomic processes are solved for k-valent ionized impurity ions.
1 ∂ ′
1 ∂
(𝑉 𝑛𝑘 ) = − ′ (𝛤𝑘 ) + 𝛾𝑘−1 𝑛𝑒 𝑛𝑘−1 − 𝛾𝑘 𝑛𝑒 𝑛𝑘 − 𝛼𝑘 𝑛𝑒 𝑛𝑘 + 𝛼𝑘+1 𝑛𝑒 𝑛𝑘+1 ,
(13)
′
𝑉 ∂𝑡
𝑉 ∂ρ
1 ∂ ′
1 ∂
(𝑉 𝑛1 ) = − ′ (𝛤1 ) − 𝛾1 𝑛e 𝑛1 + 𝛼2 𝑛𝑒 𝑛2 − 𝛼1 𝑛𝑒 𝑛1 + 𝑆0 (𝑘 = 1),(14)
𝑉 ′ ∂𝑡𝑉 ∂ρwhere, γk is the ionization rate, αk is the recombination rate and S0 is the impurity source. The particle flux Γ is thesame as in Eq. (3), and the turbulent diffusion coefficient is DZ = Di. The impurity temperatures are assumed to beequivalent to the bulk temperature. The γk, αk and radiation loss power of impurity are evaluated by using OPENADAS [12] database. In this simulation, the impurity source is set for the singly charged impurity ion, assuming thatthe ion density is n1 to have thermal velocity vth at the plasma edge. The time evolution of bulk plasmas and impurities
is solved by linking the TR and TI modules. The data of plasma density and temperature profiles are transferred fromTR to TI, while the data of radiation loss power and effective charge number are transferred from TI to TR.
𝐹(𝑠, α) =3.

Tungsten transport in negative magnetic shear mode

In this section, the results of negative magnetic shear mode calculations arepresented with the parameters referring to JET tokamak; the major radius R = 3.0 m,
minor radius a = 1.2 m, ellipticity κ = 1.5, magnetic field Bt = 3.0 T, and on-axis
heating power for electrons and ions 5MW. The CDBM model is used for turbulent
transport in this calculation. Figure 1 shows the time evolution of the plasma current
and heating power. The plasma current is steady to be 3 MA in the corresponding Lmode, but in the negative shear mode calculation, it is ramped up from 1 MA to 3
MA in 1 second. The additional heating of 5 MW is applied at the same time as thecurrent ramp-up. As shown in Fig. 2, the current ramp-up operation causes the OH
current to rise from the plasma periphery, resulting in a concave current density
profile. Corresponding to the current profile, a negative shear region is formed.
The radial profiles of the bulk plasma and impurity are shown in Fig. 3. In thenegative shear mode, a strong gradient region is formed in the ion temperature anddensity profiles. This is because the transport is improved by locally reduced
diffusion coefficient in the negative shear region. The formation of the negative
magnetic shear is transient, and the ITB moves toward the plasma center with current
diffusion, eventually approaching the L-mode state. The temperature increases,
because the heating power is greater at the plasma center, and the density decreases,
because the particle source is only set at the plasma periphery. When the strong ITB
is formed, an inward pinch is driven in the ITB region to increase the tungsten
density, resulting to give a step-like profile. This induces a greater amount of
impurity accumulation in the plasma core than in the L-mode case. Especially fromt = 2 to 3 sec, the tungsten density increases due to large bulk density gradient. For
t = 4 sec the bulk density decreases, while the strength of the temperature gradientis maintained, so the tungsten convection directs outward to reduce the tungsten
density.

Fig. 1: Time evolution of (a)plasma current and (b)heating power.

Fig. 2: Radial profiles of (a)current density and (b) safety
factor.

Fig. 3: Radial profiles of (a) plasma density, (b) ion temperature, (c) ion thermal diffusivity, (d)tungsten densityand (e) neoclassical pinch.

Fig. 4: Radial profiles of (a) plasma density, (b) electron and ion temperature, (c) electron and ion thermal
diffusivity, (d) tungsten density and neoclassical pinch.

Fig. 5: Radial profiles of (a) plasma density, (b) ion temperature, (c) tungsten density and (d) neoclassical pinch.

4.

Tungsten transport in high βp mode

In this section, the transport calculations for tungsten are carried out with reference to the high βp mode discharge
in JET [13]. The plasma parameters are R = 3.0 m, a = 0.9 m, κ = 1.5, Bt = 2.6 T and Ip = 2.2 MA. Neutral beam
injection heating PNB = 15 MW and ion cyclotron heating PICRF = 5 MW are applied as external current-driven heating
sources. As a turbulent model, the mBgB model is used, because it is an empirical model based on JET experiments.
Figure 4 shows the radial profiles of bulk plasma and tungsten impurity. The diffusion coefficient decreases due tothe influence of the magnetic shear and the E×B flow shear at ρ = 0.0 - 0.4, while there is a peak due to a numericalproblem. Inside ρ < 0.3, the effect of magnetic shear is strong, and with ρ = 0.3 - 0.4, the stabilizing effect of E×B
flow shear makes the ITB-foot position more outward. Only the ion temperature profile is compared with theexperimental data, which shows good agreement inside the ITB-foot. However, in the plasma periphery, the values
seem to deviate due to the difference in the heating or plasma density profiles. On the other hand, the density profile
of tungsten increases near the center of the plasma like the bulk density profile. This is because the diffusion
coefficient near the plasma center is smaller and the inward neoclassical pinch component is driven by the NBI
particle supply to increase the density gradient.
Calculations are performed for different NBI input power. As shown in Fig. 5, with increasing NBI power, theplasma density and temperature increase due to increased particle supply and heating, which make the transport
improvement region wider. For impurities, the inward neoclassical pinch becomes stronger and wider with higher
NBI power, and the tungsten density increases over a larger area, resulting in more tungsten accumulation.

5.

Conclusion

Tungsten transport calculations in the negative shear mode and high βp mode are performed using the integrated
code TASK. In the negative shear mode, the inward neoclassical pinch component is driven when the plasma density
gradient is maintained in the ITB region. Even in the high βp mode, it results in an increase in tungsten density in theplasma center region where the plasma density gradient is strong. Impurities may accumulate in the plasma core dueto the formation of ITB. In the future, impurity control by electron cyclotron heating, which has been confirmed tobe effective in experiments, will be investigated.
This work was supported by JST SPRING, Grant Number JPMJSP2136, and by the collaboration program of QST
on development research of the fusion DEMO reactor.

References[1] R. Dux, et al., J. Nucl. Mater. 313, 1150 (2003).
[2] W. Chen, et al., Phy. Lett. A 440, 128141 (2022).
[3] C. Angioni, et al., Nucl. Fus. 54, 083028 (2014).
[4] R. Neu et al., Plasma Phys. Control. Fusion 44, 811 (2002).
[5] https://bpsi.nucleng.kyoto-u.ac.jp/task/[6] M. Honda, et al., Nucl. Fusion 46, 580 (2006).
[7] W. Houlberg, et al., Phys. Plasmas 4, 3230 (1997).
[8] M. Yagi, et al., Phys. Fluids B 5, 3702 (1993).
[9] A. Fukuyama et al., Plasma Phys. Control. Fusion 37, 543 (1995).
[10] T. J. J. Tala, et al., Plasma Phys. Control. Fusion 43, 507 (2001).
[11] L. Garzotti, et al., Nucl. Fus. 43, 1829 (2003).
[12] ADAS website, http://www.adas.ac.uk/[13] P. Buratti, et al., Nucl. Fus. 41, 1809 (2001).

Particle Transport in Fusion Plasmas
M. Yagi
National Institutes for Quantum and Radiological Science and Technology,
Rokkasho Fusion Institute, Rokkasho, Aomori, Japan
1. Introduction
In fusion plasmas, the burn control is important issue to sustain the self-ignition condition. We need toprovide fuels externally across the last closed flux surface, at the same time, extract Helium ash fromthe plasma core. Therefore, to understand particle transport mechanism is a key to develop the burncontrol method. For fuel injections such as gas puffing or pellet injection, we expect the hollow
density profile will be formed in plasma peripheral region. The stabilization of trapped-particle modes
by the inverted density gradient was reported [1]. According to local gyrokinetic simulations, wefound the inward particle flux was driven by the inverted density gradient in the semi-collisional
plasmas [2]. The particle transport driven by the ion-mixing mode [3] was also investigated. We
analytically found that the growth rate of electron drift wave driven by the parallel thermal
conductivity (semi-collisional effect) and negative he is larger than the ion mixing mode with hi < 0limit where hi,e implies the ratio of the density gradient scale length to the ion or electron temperature
gradient scale length, respectively. However, the question is naturally arisen whether such follow
profile is sustained or not under plasma turbulent state, in other words, how long does the source ofinward particle flux survive? To answer this question, the global simulation is inevitable. In thepresent study, the global fluid simulation based on the 6-field drift reduced MHD model is performed
which is the extension of 4-field model [4].

2. 6-field drift reduced MHD model
The 6-field drift reduced MHD model is composed of the vorticity equation, the generalized
Ohm’s law, continuity equation, the parallel momentum equation, the electron temperature
evolution equation and ion temperature evolution equation. They are given by
!
!"
!
!"
*+
*"

∇$# 𝐹 + [𝐹, ∇$# 𝐹 ] − 𝛿∇# ∗ [𝑝% , ∇# 𝐹 ] = −∇∥ 𝐽 + 𝜇#% ∇'# 𝐹
(

(

(1)
(

.𝐴 − 𝛿 $ (! 𝐽0 = −∇∥ (𝜙 − 𝛿𝑝) ) + 𝛿 $ (! [𝜙, 𝐽] + 𝜂∥ 𝐽 − 4𝜇#) 𝛿 $ (! ∇$# 𝐽
"

*,

"

+ 𝛽 *" = 𝛽[Ω, 𝜙 − 𝛿𝑝) ] − 𝛽∇∥ (𝑣 + 𝛿𝐽) + 𝜂# 𝛽∇$# 𝑝

"

(2)
(3)

*-

(/)

= −∇∥ 𝑝 + [𝐴, 𝑝] + 4𝜇#% ∇$# 𝑣

*"
*1!

2 *+

*1"

2 *+

(4)
3

− 2! *" = −𝛼 1 𝛿𝛽) ∇∥ 𝐽 + 𝜖 $ 𝜒∥) ∇$∥ 𝑇) + 𝜒#) ∇$# 𝑇) − $ 𝛿𝛽) [Ω, 𝑇) ]
*"

(5)

3

− 2" *" = 𝜖 $ 𝜒∥% ∇$∥ 𝑇% + 𝜒#% ∇$# 𝑇% + $ 𝛿𝛽% [Ω, 𝑇% ]
*"

(6)
4

'7,

'7,

where 𝐹 = 𝜙 + 𝛿𝑝% , 𝐽 = ∇$# 𝐴, Ω = 2𝑟𝑐𝑜𝑠𝜃, 𝛿 = 5 6 , 𝛽) = 8# ! , 𝛽% = 8# " , 𝛽 = 𝛽) + 𝛽% ,
"

*

!

*"

(/)

= !" + [𝜙, ] represents the convective derivative and ∇∥ = ∇∥ − [𝐴, ].

The normalization is performed based on the minor radius 𝑎 and the poloidal Alfven time
6

where 𝜖 = 𝑎/𝑅/ is the inverse aspect ratio, 𝑅/ is the major radius. The energy

9-$conservation law is written as*〈ℎ〉 = −〈𝑠〉

*"

(7)

with
:

:

:

(

:

:

: +#

: ;1 #: ;1 #!

"ℎ = $ |∇# 𝐹|$ + $ |∇# 𝐴|$ + $ 𝛿 $ (! 𝐽$ + $ 𝑣 $ + $ 𝑝$ + $ 2 + $ $2! + $ $2"
"

(8)

and
(

𝑠 = 𝜇#% |∇$# 𝐹|$ + 𝜂∥ 𝐽$ + 4𝜇#) 𝛿 $ (! |∇# 𝐽|$ + 𝜂# 𝛽|∇# 𝑝|$ + 4𝜇#% |∇# 𝑣|$ +
"

<&!
2!

|∇# 𝑇) |$ +9 # <∥"
2"$<

J∇∥ 𝑇% J + 2&" |∇# 𝑇% |$"9 # <∥!
2!

$

J∇∥ 𝑇) J +(9)where 〈⋯ 〉 represents the volume integral. It should be noted that 𝑠 is given by the positive
definite form [5].
The density source and the sink term were introduced into the continuity equation (3) and theelectron temperature evolution equation (5), respectively. Here, we assume that the density
source has a sphere shape:=#?#

𝑆+ = 𝑆/ exp .− $>#0 exp .− $># 9#0(10)

where 𝜉 = Q𝑟 $ + 𝑟@$ − 2𝑟𝑟@ cos 𝜃 . For the sink term, we simply assume the adiabatic
expansion of the source, namely, 𝑆1! = −𝑆+ .
3. Simulation results

For the simulation parameters of R6F, we use 𝑆/ = 0.2, 𝑟@ = 0.8, ∆= 0.1, 𝛽% = 𝛽) =4 × 10A; , 𝛿 = 110A$ , 𝜖 = 1/3. The mesh sizes are chosen as 513 × 256 × 128. The initial
profiles are shown in Fig.1.

Fig.1

Initialtemperature

𝒏

,

and

𝒒

density𝑻𝒊 = 𝑻𝒆profiles. ITG mode is unstable
for these profiles.

10-3

t=4000edens( 0/ 0)
edens( 2/ -1)
edens( 3/ -1)
edens( 1/ 0)
edens( 4/ -1)
edens( 1/ -1)
edens( 5/ -2)
edens( 4/ -2)
edens( 11/ -3)
edens( 10/ -3)

10-5

10-7

t=4500En
10-9t=5000
10-11

10-13

0

2 103

4 103

6 103

8 103

1 104

t

Fig.2 Time evolution of fluctuating density
energy of each Fourier mode.

Fig.3 Contour plots of fluctuating
ion temperature at z=0 plane.

Fig.2 shows the time evolution of fluctuating density energy of each Fourier mode which is
given by the sixth term in Right Hand Side of eq. (8). At the early phase, ITG modes linearlygrow up and the saturation is obtained at t=4300 due to the ion temperature profile relaxation.
Then, the density source is applied at t=6400. Fig.3 shows the contour plots of fluctuating ion
temperature on z=0 plane at t=4000, t=4500 and t=5000, respectively. It is clearly observed
ITG mode structure at t=4000 and then the ion temperature relaxation occurs at t=5000. Fig.4
shows the contour plots of helical flux function, fluctuating density, total density andtemperature at t=10000.

𝜓∗𝑛" + 𝑛$𝑛$

𝑇#" + 𝑇'#Fig.4 Contour plots of helical flux function, fluctuating density, total density andtotal temperature at t=10000.

It is seen that density source does not propagate homogeneously. This is because (2,1)
magnetic island forms at r=0.5 which plays a role of barrier for density source propagation. It
is clearly observed in the flux averaged total density and temperature shown in fig.5.

2.5 10-2

1.2 10-2
t=6400t=8000t=10000-2

2 10

t=6400t=8000t=10000

1 10-2
8 10-3

1.5 10-2
Te 6 10-3n
1 10-2

4 10-3
5 10-30

2 10-3

0

0.2

0.4

0.6

0.8

1

r

0

0

0.2

0.4

0.6

0.8

1

r

Fig.5 Flux averaged total density and temperature at t=6400(dashed),
t=8000(blue), and t=10000(red), respectively.

The discontinuity point of radial derivative of the flux-averaged density appears around r=0.5
which indicates the propagation characteristics of density source is changed at this point. In
addition, we observe the density hump on the mid-plane in Fig.4, however, it does not appear
in the flux-averaged total density profile at t=10000.
4. Summary and discussion
We have developed 6-field drift reduced MHD model to investigate global particle transportwith density source. The local density on the mid-plane shows density hump, however, forflux-averaged density profile, the hump does not form at the present time scale. We speculatethat it is too short to equilibrate the density profile. In addition, well-developed turbulence(electron drift wave etc.) can not be captured for the present simulation. More long-time
simulation run might be necessary. It is left for future work.
Acknowledgements
The author would like to thank Drs. N. Miyato, A. Matsuyama, N. Kasuya and T. Takizuka for
useful comments. The computations were carried out on JFRS-1 supercomputer at
IFERC-CSC and on SGI HPE8600 supercomputer at QST and JAEA.
References[1] W. M. Tang, et. al., Phys. Rev. Lett. 35 (1975) 660-663.
[2] M. Yagi, et. al., ‘Revisit of ion-mixing mode transport’, 14th Japan-Korea Workshop on
Modeling and Simulation of Magnetic Fusion Plasmas, 2019 June 27-28, RIAM Kyushu
University.
[3] B. Coppi and C. Spight, Phys. Rev. Lett. 41 (1978) 551-554.
[4] R. Hazeltine and J. Meiss, Phys. Rep. 121 (1985) 1-164.
[5] M. Yagi and W. Horton, Phys. Plasmas 1 (1994) 2135-2139.

Study of turbulent transport in trapped electron and ion temperature
gradient modes in high collision regime for tokamak plasmas
S. Toda, M. Nunami, and N. Kasuya1
National Institute for Fusion Science, Oroshi-cho 322-6, Toki, Gifu, 509-5292, Japan
1 Research Institute for Fusion Science, Kyushu University, Kasuga-koen 6-1, Kasuga, Fukuoka, Japan

The purpose of this study is to quantify turbulent transport in toroidal plasmas. Microscopic
plasma instabilities are studied using electromagnetic gyrokinetic simulations. Linear calculations
are used to investigate the type of instability and to evaluate the conditions for instability excitation. The saturation level of the instability is determined from nonlinear calculations to obtain the
turbulent transport value. The plasma experiments at the PLATO device at the Research Institute for Applied Mechanics, Kyushu University [1] have been started. The plasma proﬁles in the
PLATO are predicted using the integrated code, the TASK code [2]. The gyrokinetic analysis havebeen performed for the plasma proﬁle and magnetic ﬁeld conﬁguration which were predicted by
the TASK, as an example in the PLATO. The instabilities are shown to be exited by the trappedelectron mode (TEM) because of the increase of the density gradient [3] where the characteristicfrequency is between the ion and electron bounce frequencies [4, 5] and ion temperature gradient(ITG) mode in the predicted parameter regions for the PLATO. Previously, an analytic model,
the ﬂuid-based model and the gyroﬂuid model for the ITG-TEM analysis have been reported.
The local gyrokinetic simulation has been performed for the ITG-TEM plasmas. Although theglobal gyrokinetic codes have been developed and performed, a local gyrokinetic simulation using
the model collision operators [6, 7] is performed in this study. First, the type of the instability
and the condition under which the instability occurs are studied by the linear gyrokinetic simulations. Next, the nonlinear simulation demonstrates the time evolution of the electrostatic potential
ﬂuctuations. The results of the simulations obtained by using the Sugama collision operator [6]are compared with those by the Lenard-Bernstein collision operator [7]. In the previous studies,
the reduced models are constructed in the adiabatic electron condition and the kinetic electron
condition [8] for the ion temperature gradient (ITG) turbulence in helical plasmas. These reduced
models, which are constructed by the linear simulation results, quickly reproduce the nonlinear
simulation results for the turbulent transport, thus, these reduced models enable us to simulate
the integrated (dynamical) transport simulation, when the turbulent transport is evaluated by the
reduced models in the dynamical transport simulation [9]. The simulation study for the TEM and
ITG mode plasmas in this article is the basis when the reduced models for the TEM and ITG
mode turbulence in tokamak plasmas will be constructed.
The turbulence driven by microinstabilities for the predicted plasmas in the PLATO is studied.
The plasma experiments at the PLATO device [1] have been started. The plasma proﬁles [2]predicted by the TASK code in the PLATO are used. Toroidal plasmas with concentric circular
ﬂux surfaces are studied. For simplicity, the proﬁle of the safety factor, q(ρ) = 1 + 2ρ2 is set.
The detailed analysis for the safety factor in the PLATO has been performed [10]. In this study,

the GKV code [11] is used for the local ﬂux tube gyrokinetic simulation. The normalized collision frequencies, νe∗ and νi∗ are shown as 0.22 and 8.8 at ρ = 0.65, respectively. Therefore, theelectrons form the banana orbit and the ions do not. The normalized electron and ion radii areabout ρ∗e = 3.9 × 10−4 and ρ∗i = 6.5 × 10−2 at ρ = 0.65, respectively. Therefore, the GKV code isperformed in the unit of R/LTe , R/vte and ky ρe .
The linear gyrokinetic simulation is performed in tokamak plasmas. The linear simulation results
using the Sugama (S) collision operator and the Lenard-Bernstein (LB) collision operator arecompared. The dependences of the linear growth rate, γ on the normalized poloidal wavenumber, ky ρe are studied at ρ = 0.47, 0.65 and 0.81. The real frequency, ω as the function of ky ρe
are also investigated. The instabilities predicted in this study are driven by the trapped electronmode (TEM) and the ion temperature gradient (ITG) mode at ρ = 0.47 and ρ = 0.65. The real
frequency changes from the electron to ion diamagnetic drift motion direction at ρ = 0.47 andρ = 0.65 when the poloidal wavenumber increases. The TEM is unstable, where the normalized
electron plasma frequency is less than the unity [4, 5]. The phenomena such as the transition fromthe TEM to the ITG mode are obtained in the dependence of the real frequency on the poloidal
wavenumber, specially at ρ = 0.65 using the S collision operator. The TEM mode is predicted
to unstable at ρ = 0.81. Using the S collision operator, the ITG mode is unstable in the region
0.026 ≤ ky ρe < 0.100 at ρ = 0.47 and in the region 0.030 ≤ ky ρe < 0.100 at ρ = 0.65. On theother hand, the ITG mode is stable in these poloidal wavenumber regions, when the LB collision
operator is used.
Nonlinear gyrokinetic analysis (results are ⟨studied for⟩ the) time evolutions of the squared turbulent potential ﬂuctuation, T = Σk̃x ,k̃y ̸=0 |ϕ̃k̃x ,k̃y |2 /2 and the squared zonal ﬂow potential,⟩ )(⟨Z = Σk̃x ,k̃y =0 |ϕ̃k̃x ,k̃y |2 /2 at ρ = 0.47, ρ = 0.65 and ρ = 0.81, where ϕ̃ is the electrostatic
potential ﬂuctuation. The electron and ion heat diﬀusivities are obtained as 29m 2 /s and 32m2 /s
at ρ = 0.65, respectively. The particle diﬀusivity is shown as 19m 2 /s. The nonlinear simulation
using the S collision operator will be performed.
This work is partly supported by the NIFS Collaboration Research Programs, NIFS22KIST022(Plasma Simulator), the JSPS KAKENHI Grant No. 21K03514 and the Collaborative Research
Program of Research Institute for Applied Mechanics, Kyushu University, 2022FP-02.

[1] A. Fujisawa, AIP Conference Proceedings 1993, 020211 (2018)[2] N. Kasuya, The 16th Burning Plasma Simulation Initiative (BPSI) Meeting (2018)[3] B. Coppi and G. Rewoldt, Phys. Lett. 49A, 36 (1974)[4] B. B. Kadomtsev and O. P. Pogutse, Reviews of Plasma Physics, edited by M. Leontovich(consultants Bureau, New York, 1970), Vol. 5[5] W. M. Tang et al., Phys. Rev. Lett. 35, 660 (1975)[6] H. Sugama, T. -H. Watanabe and M. Nunami, Phys. Plasmas 16, 112503 (2009)[7] A. Lenard and I. B. Bernstein, Phys. Rev. 112, 1456 (1958)[8] S. Toda et al., Phys. Plasmas 26, 012510 (2019)[9] S. Toda, M. Nunami and H. Sugama, Plasma Phys. Control. Fusion 64, 085001 (2022)[10] S. Tomimatsu et al., Plasma and Fusion Research 15, 1403052 (2020)[11] T. -H. Watanabe and H. Sugama, Nucl. Fusion 46, 24 (2006)

Several Low-Cost Techniques for Global Plasma Simulation with
Particle-in-Cell Model
T. Takizuka 1, K. Ibano 1, S. Togo 2
1
2

Graduate School of Engineering, Osaka University, Suita 565-0871, Japan
Plasma Research Center, University of Tsukuba, Tsukuba 305-8577, Japan

Kinetic particle model is valuable even for global simulations of low-collisionality fusion plasmas,
although the computational cost is very high. We have been proposed several techniques tolower the computational cost for global plasma simulations with Particle-in-Cell model;“Weighted Particle model” for hugely multiplying particle number, “Double Leap-Frog
method” for larger time step, “Ingenious plasma simulation model” for longer spatial grid size,
etc. Even when the above models are used, the numerical heating is unavoidable especially forsmaller particle number per cell. Here we propose a “Numerical Heating Correction method”
and “Increasing/Decreasing Particle Number techniques”.

1. Introduction
Numerical simulations have been widely carried out for various fields of plasma science. Global
plasma simulations, exceeding local simulations, are important especially for fusion plasma R&D.
Since kinetic effects are essential in low-collisionality plasmas, the kinetic particle model is valuable
even for global simulation [1]. The computational cost of the particle simulation, however, is veryhigh. For example, a three-dimensional (3D) particle simulation of a tokamak plasma (volume Vp ~10 m3 and time scale Tsim ~ 10 ms) including the self-consistent electron dynamics needs large amountsof super-particle number Ntot ~ 1013 and time steps Kstep ~ 108 (grid size Δ ~ 1 mm, particle number percell Ncell ~ 103, and time step Δt ~ 10–10 s). Resultant calculation steps become very huge Kcal ~ 1024
(for operation per step ~ 103). Using 400 PFLOPS “Fugaku” system, a single run may require ~ 1 year.
To reduce the computational cost, the super-particle concept has been applied in the Particle-in-Cell(PIC) modeling. In the above tokamak plasma with electron density ne ~ 1019-20/m3, the real particle
number Nreal ~ 1020-21 is reduced to the super-particle number Ntot ~ 1013. Since Coulomb collisionnature is deformed in the PIC model, a “binary collision model” based on a Monte-Carlo techniquewith low cost was developed [2], and has been applied widely to the global plasma simulations.
When we study with the particle simulation the vapor shielding effect at the ELM or disruption
event, large amount of impurities are emitted from the wall and the particle number is drasticallyincreased. In order to suppress the corresponding increase of the super-particle number, a “weighted
particle model” was introduced to the PIXY code and the computational cost was much reduced [3].
Leap-frog (LF) method is simple to follow the gyration motion of a charged particle with energy andmagnetic moment conservation and accurate drift. The simulated Larmor radius, however, becomesvery large for ΩΔt >> 1 (Ω is cyclotron frequency). If the LF method is applicable to large ΩΔt case,new particle simulation codes can be more easily developed, such as the PIXY of multi-species(electron, H+/D+/T+, impurities) code, in which the same subroutine is utilized for all particle motions.
We proposed therefore a “double leap-frog method” [4]. By repeating the LF steps doubly andadopting the averaged velocity to advance the particle position, the Larmor radius is kept real
independently of the ΩΔt value.
When the size of a super-particle (or the grid size Δ) becomes much larger than the Debye length λD,
the explosive heating called “finite grid instability” occurs with a growth rate γ of the order of plasma
frequency ωp, and decent simulation results are never obtained [5]. In a tokamak fusion plasma, λD isvery short ~ 0.1–0.01 mm for ne ~ 1019–20 /m3 and electron temperature Te ~ 102–4 eV. This length ismuch shorter than a possible grid size Δ ~ 1 mm in a global simulation. We therefore proposed an“ingenious plasma simulation model” stably carried out even for the longer Δ [6,7]. An artificial

permittivity ε* much greater than ε0 (permittivity of free space) is introduced to the Poisson equation,∇(ε*∇φ) = −ρ (φ is electrostatic potential and ρ is charge density), and λD* is artificially elongated by(ε*/ε0)1/2 times. A large-scale system can be treated with a reasonable number of spatial grids. A
dielectric tensor is adopted for the 3D global simulation of tokamak plasmas in the cylindricalcoordinates (R, θ, Z). A diagonal element εθθ for toroidal axis θ is set much larger than those for thepoloidal plane, εRR = εZZ = ε*. Number of toroidal meshes Nθ can be reduced enough smaller than thatof poloidal-plane meshes, NR and NZ, in comparison with Nθ > NR, NZ for the usual case of isotropic ε*.
Even for Δ ~ λD* of the “ingenious model”, the numerical heating is unavoidable in the PIC
modeling. This heating rate is in inverse proportion to the super-particle number per cell, Ncell. In spiteof the large numerical heating risk, it is convenient to save Ncell ~ 100 for the reasonable computation
cost of multi-dimensional particle simulation. The cost of the former 3D tokamak simulation (Vp ~10 m3 and Tsim ~ 10 ms) can be reduced much by the following numerical condition; ΔR,Z ~ 2 mm, Δθ ~10 mm, Δt ~ 10–9 s, Ncell ~ 100, and Ntot ~ 1010. Resultant run of Kcal ~ 1020 requires ~ 1 day by 10
PFLOPS “Raijin” system. In the next section, we briefly review the numerical heating, and propose a“numerical heating correction method”. In section 3, we also propose “increasing/decreasing
particle number techniques”.
2. Numerical heating correction method [8]In the PIC model, the electric field E is recorded at the grid points, and nonphysical aliases of E
arise on an off-grid-point particle. This aliasing effect causes the finite-grid instability for Δ >> λD (orλD*). Even for Δ ~ λD*, the aliasing E noise suffers the electron velocity parallel to the magnetic fieldB, and causes the numerical heating in the parallel component of electron energy. Major source of thenumerical-heating field is of the short wave-length statistical fluctuation, whose amplitude is in inverse
proportion to (Ncell)1/2; e δE Δ /Te ~ fA (Δ/λD*)2/(Ncell)1/2. Aliasing factor fA can be ~0.1 for the linear-shape
PIC. Electron velocity is diffused by this noise field, and the diffusion coefficient in the velocity space
is Dv ~ (e δE /me)2 Δteff , where Δteff ~ {4(Δ /vte)2 + Δt 2}1/2 is the decorrelation time (vte is electron thermal
speed). The numerical heating rate is then evaluated as γN (≡ d lnTe /dt) ~ Dv /vte2 ~ ωp* fA2 (Δ/λD*)3/Ncell
for Δt < Δ /vte [1]. Simulation result by Ueda et al. was γN ≈ (1/1600) ωp (Δ/λD)3/Ncell in a periodic 1D
system [9]. This unphysical increase of Te becomes quite large for a confined plasma, and interrupts along-time scale global simulation. The higher-order shaping of a super-particle (quadratic spline andcubic spline) can reduce the heating rate [10], but still cannot fully suppress the numerical heating.
Because the linear spline (PIC) is simple for the multi-dimension system and is preferable for theboundary treatment in the open system, we do not adopt the higher-order shaping while we consider acorrection method of the numerical heating in the PIC simulation. Complete suppression of thenumerical heating is accomplished by the following correction procedures.
Correction procedures for every time step
( i ) Simply smoothed electric field parallel to B, E|| = (E⋅B)B/B2, is calculated at a grid point j;〈E||j 〉 = (E||j+1 + 2E||j + E||j–1)/4(1)

Fluctuation energy of “smoothed” field Σj 〈E||j 〉2 becomes small by ~1/3 compared with the “original”
fluctuation energy Σj (E||j )2.
(ii) Each electron velocity v(t) is advanced to v*(t+Δt) with the “original” E. The perpendicular component is determined, v⊥(t+Δt) = v*⊥(t+Δt), at this step.
(iii) Each electron velocity v||(t) parallel to B is advanced to v#||(t+Δt) with the “smoothed” 〈E||〉;v#(t+Δt) = v*(t+Δt) – (e/me) (〈E||〉 – E||) Δtv#||(t+Δt) = {v#(t+Δt)⋅B}B/B

2

(2)
(3)

(iv) Total energy increments (including ion contribution), ΔW* and ΔW#, are calculated for the“original” E step (ii) and for the “smoothed” 〈E||〉 step (iii), respectively.
( v ) Input and output powers in the total system (hot particle supply, ohmic heating, external heating,fusion power, convection/conduction loss to the wall, radiation loss, atomic & molecular reaction loss etc.)
are counted. Physical heating power Pin should be accurately evaluated. Physical energy balance is thenknown; ΔW0 = (Pin – Ploss) Δt. After that numerical heating energies, ΔWN* and ΔWN# , are calculated;

ΔWN*,# = ΔW*,# – ΔW0(ΔWN* > ΔWN#)(4)(vi) Assuming the numerical heating comes from only electrons, the parallel electron energy is cooledby keeping the direction of the parallel velocity asv||2 = v*||2 – (v*||2 – v#||2) ΔWN*/(ΔWN* – ΔWN#)(5)

v||(t+Δt) = (v||2/v#||2)1/2 v#||(t+Δt)(6)

Total excess energy is fully removed in the entire system. Note that a simple cooling w/o the“smoothed” 〈E||〉 procedures causes the numerical heat transport.
(vii) Spatial position is moved (LF method); r(t + 3Δt/2) = r(t + Δt/2) + v(t+Δt) Δt .
(viii) Collisional change is added by using a binary collision model (Takizuka-Abe [1]).
Complete suppression of the numerical heating in the entire system can be accomplished by the above
correction method.
3. Increasing/decreasing particle number techniques
In order to perform a global plasma simulation with a convenient number of super-particles, it isefficient to vary arbitrarily the particle number during a simulation run. A long-duration simulation
could be carried out with Ncell ~ 102, and an accurate result could be obtained in a short period with Ncell ~103. Further reduction of Ncell ~ 10 is desirable in the pre-simulation to obtain an initial equilibrium state.
Previous PARASOL 2D simulation, the particle number was multiplied by 10 times at once from thepre-simulation (w/o self-consistent E) to the full simulation [11]. Now we propose “increasing/decreasing particle number techniques”, where we employ the “doubling/halving number processes”.
Figure 1 shows the temporal evolution of the total particle number in the system, Ntot. During theincreasing/decreasing phase, Ntot is doubled/halved at a time step. A short plateau (~102 Δt) is put after
a doubling/halving process to suppress the unexpected fluctuation induced by this process. It takes
seven doubling/halving processes (27 = 128) and ~103 Δt duration for the increasing/decreasing of Ntot
by hundred-fold/one-hundredth. A certain simulation result is obtained in a long plateau (>103 Δt).
During the plateaus, physical increase/decrease of Ntot is naturally allowed.

Fig. 1 Temporal evolution of total
particle number in the system, Ntot.
Series of doubling/halving processes,
and short/long plateaus are seen.

“Doubling process”
A super-particle with (r0, v0, w0) is doubled to two super-particles with (r1, v1, w1) and (r2, v2, w2),
respectively. The weight of doubled particle is half of the first particle, w1 = w2 = w0/2. Their new
velocities and positions are given as
v1 = v0 + (v0 + v1) × Θ p , r1 = r0 + v1Δtp(7)

v2 = v0 – (v0 + v2) × Θ p , r2 = r0 – v0Δtp(8)

where Θ p = ΩΔtp /2 and (B⋅v0)Δtp /B ~ Δ. Here the forward LF method, Eq. (7), and backward LF
method, Eq. (8), are utilized to separate the spatial positions of doubled particles. A binary collision
between v1 and v2 takes place with a scattering angle ~ π/4. The energy is conserved before and after adoubling process. Particle diffusion along B is naturally described here to reduce the source ofnumerical-heating E of the short wave-length statistical fluctuation.
“Halving process”
Two super-particles with (ra, va, wa) and (rb, vb, wb), respectively, nearby are chosen as a pair. One
remaining particle is selected, a or b (hereafter “a” for simplicity). The position ra and parallel velocity
va|| are taken over by the new super-particle, r1 = ra and v1|| = va||. The weight of the new particle is w1= wa + wb. Although the momentum is not conserved, w1v1 ≠ wava + wbvb, the energy conservation is
maintained by modifying the perpendicular velocity v1⊥ asv1⊥2 = (wava2 + wbvb2)/(wa + wb) – v1||2
2

2 1/2

v1⊥ = (v1⊥ /va⊥ )

va⊥

(9)
(10)

The momentum conservation is also possible, if necessary, by applying the cell-average correction in
the “weighted-particle collision model” [12]. Halving pair particles are chosen mainly in the same
cell, but the inter-cell pair is allowed. Last one particle not paired in the system if exist is 50%-life (w1= 2wa) or 50%-death (w1 = 0).
4. Summary
Kinetic particle model is valuable even for global simulations of low-collisionality fusion plasmas,
although the computational cost is very high. We have been proposed several techniques to lower thecomputational cost for global plasma simulations with Particle-in-Cell model; “weighted particle
model” for hugely multiplying particle number [3], “double leap-frog method” for larger time step [4],“ingenious plasma simulation model” for longer spatial grid size [6,7], etc. Even when the above
models are used, the numerical heating is unavoidable especially for smaller particle number. A“numerical heating correction method” is proposed for longer time-scale global plasma simulation
with smaller particle number [8]. “Increasing/decreasing particle number techniques” are also
proposed for the efficient simulation with arbitrary particle number. These low-cost techniques will beapplied and tested during our development of the PIXY code (multi-dimensional PIXY for the mirror
confinement system) in the near future.
References[1] T. Takizuka, Plasma Phys. Control. Fusion 59, 034008 (2017).
[2] T. Takizuka, H. Abe, J. Comput. Phys. 25, 205 (1977).
[3] K. Ibano, Y. Kikuchi, S. Togo, Y. Ueda, T. Takizuka, Nucl. Fusion 59, 076001 (2019).
[4] T. Takizuka, K. Ibano, S. Togo, Plasma Fusion Res. 16, 1203100 (2021).
[5] R.W. Hockney, J.W. Eastwood, Computer Simulation Using Particles (New York, McGraw-Hill, 1981).
[6] T. Takizuka, K. Ibano, M. Yagi, Plasma Fusion Res. 13, 1203088 (2018).
[7] T. Takizuka, K. Ibano, M. Yagi, Plasma Fusion Res. 14, 1203091 (2019).
[8] T. Takizuka, K. Ibano, S. Togo, “Correction of Numerical Heating in Global Plasma Simulation with
Particle-in-Cell Model”, presented at 27th Int. Conf. Numerical Simulation of Plasmas (Nagoya;
Online, 2022), P-6; https://www.p.phys.nagoya-u.ac.jp/27thicnsp/program.html .
[9] H. Ueda, Y. Omura, H. Matsumoto, T. Okuzawa, Comput. Phys. Commun. 79, 249 (1994).
[10] H. Abe, N. Sakairi, R. Itatani, J. Comput. Phys. 63, 247 (1986).
[11] T. Takizuka, K. Shimizu, N. Hayashi, M. Hosokawa, M. Yagi, Nucl. Fusion 49, 075038 (2009).
[12] A. Tanaka, K. Ibano, T. Takizuka, Y. Ueda, Contrib. Plasma Phys. 58, 451 (2018).

磁化プラズマにおける巨視的構造診断のためのシミュレーション
Global structure and its diagnostic simulation in magnetized plasmas
N. Kasuya
RIAM, Kyushu University, Kasuga, Fukuoka 816-8580, Japan
E-mail: kasuya@riam.kyushu-u.ac.jp
1. Introduction
Plasma turbulence has large effects on particle and heat transport in magnetically-confined
fusion plasmas, and it is important to clarify its 3-D feature in real magnetic geometries intorus experiments, because various global and meso-scale structures are formed. Synthetic
diagnostics use simulation data to represent experimental measurements by simulating
measurement signals [1]. Numerical simulations can show global features of mean profile
modifications and fluctuation excitations, and quantitative comparison between experiments
and simulations is possible. Several kinds of simulations use different models, depending on
the simulation targets, so comparison between models is also important to capture various
characteristics in magnetized plasmas. In addition, combination of multiple simulation results
can give complex plasma characteristics including several physical mechanisms. We havecarried out some numerical diagnostics [2,3], and then newly develop these routines into an
integrated numerical diagnostic platform for turbulence and structural formation analyses,
which can be applied to multiple simulation data that have different data formats [4]. The
concept of the integrated diagnostic simulator is explained and is applied to observations of aformed mean profile and fluctuation distribution in real magnetic geometries of fusion devices
here. Effectiveness of this framework is shown for detecting global features by analyzing
several kinds of numerical simulation data.
2. Integrated turbulence diagnostic simulator
The integrated turbulence diagnostic simulator (iTDS) is a framework for making global
analyses in real magnetic geometries of fusion devices [4]. iTDS connects global simulations
and experimental observations. Various analysis routines for several experimental diagnostics
are developed on this platform. The analysis consists of the following processes;i) 3-D simulation of plasma instabilities and transportii) Plasma device parameter set and adjustment of the data format (if necessary)iii) Simulation to obtain signals by experimental diagnostics
In process i), several codes can be used to obtain 3-D plasma profiles and time evolution offluctuations. Global codes are used for calculations in torus plasma geometries. In general,
particular models do not always include all the necessary quantities, so lacked information is
supplemented with other calculation or assumption, if necessary in process ii). In addition,
spatial coordinates in simulation codes may be different with ones used for numericalmeasurement, so data transfer is necessary in that case. Data mapping to the common format
is often used, but here we introduce the object-oriented method in accordance with each data
format to gain the maximum information. The data transfer is limited to the minimum. Thenthe numerical measurements are carried out in process iii). The basic component to capture
3-D dynamics from simulation data is extracting the local value at the arbitrary position fromdiscreet data by interpolation. The combination of this basic routine can give multi-point
observation along a line of the sight of each diagnostic. Here the focus is mainly set on ageometric effect from magnetic configurations. Applications of iTDS is explained in thefollowings to show mean profile and fluctuation measurements in torus experimental devices.

3. Steady profile measurement
Heavy ion beam probe (HIBP) is powerful diagnostic to give information of several kindsof quantities [5]. Single trajectory of the detected ion can be selected with given injection anddetection conditions, which determines the observation position. Plasma fluctuations andelectromagnetic fields inside the plasma can be measured simultaneously without disturbing
the plasma. Calculation of 3-D trajectory of heavy ions, including non-circular plasma shapeand experimental arrangement of diagnostics, is carried out by synthetic HIBP (sHIBP) in
iTDS for LHD [6]. Mean radial 0 and 3-D perturbation 1 profiles of the potential are
included in the 3-D magnetic configurations. The 3-D variation of the potential is calculated
by global drift kinetic code FORTEC-3D, which is given by the neoclassical process [7]. This
3-D variation is predicted to affect the impurity neoclassical transport. The observation points
depend on the injection conditions of the heavy ions, so the profile scan can be made. The
deviation of the potential values from the radial mean comes from 1, which depends on theobservation point. For assessment of observability in the experiment, the intensity of thedetected beam is important, which is also evaluated in sHIBP. Here beam attenuation is
calculated by using ionization cross sections from Lotz’s empirical formula [8]. The detectedbeam intensity is predicted to be larger than the noise level given by the detector.
4. Fluctuation measurement
Phase Contrast Imaging (PCI) utilizes the laser scattering by density fluctuations [9].
Interference between scatted beam and non-scatter beam gives the fluctuation image on theimaging plate. The line integration is obtained as the observed signal. Local information isneeded to be deduced from the line integration, and in helical plasma variation of themagnetic field direction can be used for the separation, because the fluctuation is elongatedalong the field line at each position [10]. Synthetic PCI (sPCI) is carried out to make thisseparation process of the spectrum using 3-D fluctuation data in iTDS. A fluctuation pattern is
obtained by gyrokinetic simulation in LHD [11], and the magnetic shear method is applied forcomparison with LHD experimental result [3]. Turbulence analyses using a torus global
simulation are also carried out.
5. Sumarry
Synthetic diagnostics, which give simulations of experimental signals by using global
simulation data of 3-D plasma turbulence, are carried out for precise comparison between
experiments and simulations. iTDS has been developed as the platform for synthetic
diagnostics by including 3-D magnetic geometry. Several kinds of diagnostics as HIBP and
PCI provide measurements for wide range of plasmas. Diagnostics using plasma waves, andmeasurements of plasma distribution function are also useful, which are future tasks.
Acknowledgements
This work was partially supported by JSPS KAKENHI JP17H06089, JP20K03905, and thecollaboration program of NIFS (NIFS22KIST029) and of RIAM of Kyushu University. Some
calculations were carried out using the JFRS-1 supercomputer system at IFERC-CSC in QST
and Plasma Simulator supercomputer system in NIFS.
References[1] C. Holland, et al., Phys. Plasmas 16 (2009) 052301.

[2] N. Kasuya, et al., Plasma Sci. Tech. 13 (2011) 326.
[3] N. Kasuya, et al., Nucl. Fusion 58 (2018) 106033.
[4] N. Kasuya, et al., Proc. BPSI meeting (2022) PA-5.
[5] T. Ido, et al., Rev. Sci. Instrum. 77 (2006) 10F523.
[6] R. Yoshihara, et al., BPSI meeting (2023) PA-5.
[7] K. Fujita, et al., Nucl. Fusion 61 (2021) 086025.
[8] W. Lotz, Astrophys. J., Suppl. Ser. 14 (1967) 207.
[9] H. Weisen, Rev. Sci. Instrum. 59 (1988) 1544.
[10] K. Tanaka, et al., Rev. Sci. Instrum. 79 (2008) 10E702.
[11] M. Nunami, et al., Phys. Plasmas 19 (2012) 042504.

Carbon impurity modeling with the superstage model in
TASK/TX
M. Honda, A. Matsuyama* and Y. Homma*
Kyoto University, Nishikyo-ku, Kyoto 615-8530 Japan*National Institutes for Quantum Science and Technology, Rokkasho, Aomori 039-3212 JapanTASK/TX features and its extension with impurities
We are developing a fluid-type transport code TASK/TX [1,2] for simulating theevolution of physical quantities in torus plasmas. Unlike the diffusion-type equations,which are the governing equations of conventional transport codes, TASK/TX istransformed into an equation system that focuses on transport phenomena by taking adrift ordering for the two-fluid equation system. Furthermore, by taking a magnetic
surface average, it becomes a one-dimensional fluid equation system. Also, in order tohandle neoclassical transport self-consistently, it solves the !1, 𝒗, !"! " , !"𝑣 " 𝒗& moment
equations of the kinetic equation up to the heat flux order, which is not usually handled
in the two-fluid equation system. TASK/TX is characterized by self-consistency between
physical quantities, and thanks to its characteristics, the quasi-neutrality condition or the
non-ambipolar condition are not imposed, which are the usual assumptions inconventional transport codes.
Recently, TASK/TX has been extended to handle impurities. To verify whether itworks with impurities, we conducted a check of ambipolar transport in both the Matrix
Inversion, which has been updated with a fully linearized collision operator implemented[3], and TASK/TX. Resultantly, it was confirmed that ambipolarity was held for bothcodes and also both results were in good agreement.

Superstage modeling in TASK/TX
As a major difference from hydrogen isotopes as bulk ions, impurities can exist in
multiple charge states within the plasma. Carbon impurities are often assumed to be in
the fully ionized hexavalent state, but this assumption may not always be valid in the edge

region. Also, assuming that the walls and divertor surfaces are made of carbon, physicaland chemical sputterings caused by ions colliding with them, reflection and pumping of
impurities at the divertor surfaces, and impurity injection are possible sources and sinksthat should be considered, but in any case, impurities enter the plasma as neutral particles.
Therefore, to model impurity sources appropriately, it is necessary to consider the multistage ionization process from neutral state to fully ionized state and the charge exchange
process. So far, impurities have been able to handle only one charge state in TASK/TX,
so impurity particle source modeling was also incomplete. By introducing the super-stage
model [4], it has become possible to easily take into account the effects of the multicharge
states from neutral carbon supply to full ionization. Solid lines in figure 1 show the time
evolution of the carbon densities of each valence and dashed lines shows the time
evolution of the densities, assuming that each valence takes a quasi-equilibrium state andgrouping all charge states of carbon into one superstage carbon. In the figure, even whensuperstaged, C+6 converges to the same density as when the reaction rate equations weresolved for all charge states, after a time longer than 10-3 second. In other words, sincecarbon is ionized from neutral to full ionization state within the time step width normallyused in TASK/TX, the application of the superstage model that treats 0-5 valence carbon
as pseudo neutral carbon is shown to be valid in the transport process that TASK/TX
focuses on.

References[1] M. Honda and A. Fukuyama, J. Comput. Phys. 227, 2808 (2008).
[2] M. Honda and A. Fukuyama, Comput. Phys. Commun. 208, 117 (2016).
[3] M. Honda, Phys. Plasmas 21, 092508 (2014).
[4] A. Foster, On the Behaviour and Radiating Properties of Heavy Elements in Fusion Plasmas,
PhD thesis, University of Strathclyde (2008).

Fig. 1: At Te = 200 eV, solid lines represent the time evolution of the densities of each valence whencarbon impurities from 0 to 5 valence are superstaged, and dashed lines represent the time evolution of
the densities of each valence when the reaction rate equations are solved.

Analysis of turbulent particle fluxes by using 5-field RMHD codein torus plasmas
H. Todoroki1, N. Kasuya1,2, M. Yagi3
IGSES, Kyushu Univ.1, RIAM, Kyushu Univ.2, QST3
1. Introduction
Fueling to the core plasma is one of the key issues to maintain a high fusion reaction rate in fusion
plasmas. Pellet injection is promising as a fueling method because, in principle, it can efficiently supply
particles to the core plasma. In fueling by pellet injection, a high-density region is created near theplasma edge due to pellet ablation, resulting in a density distribution with an inversed hump-shaped
gradient in the radial direction. In the previous study [1], inversed density gradients are known to
produce turbulent particle pinch. Furthermore, it is necessary to understand the transport mechanismsthat are not only limited to local models but are also valid in global models.
In order to elucidate the mechanism of the particle pinch effect, which is effective in fueling to theplasma center by pellet injection into the fusion plasma, we perform turbulence simulations with particle
transport by introducing a distribution with an inversed density gradient near the plasma edge. Particle
convective transport is evaluated to clarify the time evolutions of the excited instabilities, when an
inversed density gradient is introduced into tokamak plasmas.
2. Numerical model
Reduced MHD code R5F [2] is used for the simulation. The R5F code is a global simulation code using
a 5-field reduced MHD model with neoclassical viscous effects. In this model, there are five variables to besolved: the potential 𝜙, the vector potential parallel to the magnetic field 𝐴, parallel velocity 𝑉, electron
density 𝑁, electron temperature 𝑇 . For these five fields, the vorticity equation, Ohm's law, the parallel ion
momentum equation, continuity equation, and the electron temperature evolution equations are used to solve
the time evolution in 3-D space [3]. Each variable is normalized by using the minor radius in the spatial scales
and the poloidal Alfven time in the temporal scales. The following tokamak plasma parameters are used forthe simulations; major radius 𝑅
skin depth 𝛿

1.5 𝑚, minor radius 𝑎

𝑁

10 , normalized ion

10 . The inversed density gradient near the plasma edge due to pellet injection is given

with the density source peak (density intensity 𝑆𝑁

0.5 𝑚, plasma beta 𝛽𝑆

exp) initially at the radial position 𝑟

𝑟 as

.

3. Linear phase results
In this study, we analyzed the cases with the initial density profiles for source density intensities of 𝑆
1.0 and 1.5, and source locations of 𝑟

0.8. Figure 1 shows the initial density and safety factor profiles.

We call the region of positive and negative density gradients near the position 𝑟𝑟 the inversed andnormal gradient regions, respectively. As shown by the dashed lines in Fig. 1, there are resonant modes, forexample, 𝑚, 𝑛

25,8 and

36,14 in the normal and inversed gradient regions, where 𝑚 and𝑛 are the poloidal and toroidal mode numbers, respectively. Linear instabilities with ballooning structures
localized at the low field side are excited. Direction and magnitude of the particle fluxes induced by thesemodes are evaluated. Figure 2(a) shows the radial profile of particle fluxes, whose magnification in the
inversed gradient region is shown in Fig. 2(b). In Figure 2(a), the 𝑛

8 component, which is strongly

excited in the normal gradient region, is dominant overall to drive an outward particle flux. On the other hand,Fig. 2(b) shows that there is contributor to induce an inward particle flux in the inversed gradient region.

Fig.1 Initial density profiles with𝑆

Fig.2 (a) Radial profile of particle fluxes in the linear

1.0 and 1.5. The safetyphase at t = 65. The magnification of the inversed

factor profile is also shown.

gradient region is also shown in (b).

The mode excited in the negative density gradient region is the resistive electron drift wave mode, and thedependence of its particle pinch effect on the density gradient is evaluated. Figure 3 shows the results ofcalculation for the linear growth rate in the inversed gradient region (ρ = 0.76) using electron drift wave mode
theory. An instability boundary exists at the density gradient length 𝐿 ~with 𝑆

1.0 and 1.5 have 𝐿

consistent that the case with 𝑆

8 and

8. The turbulence simulations

13 in the negative gradient region, respectively, so it is

1.5 is more unstable to excite the electronic drift wave modes, though

the case with 𝑆

1.0 is less unstable. Therefore, the time evolutions of the mode are also different toeach other. The 𝑛

14 component is excited in the inversed gradient region as the electron drift wave

mode, and its particle flux increases in time in the unstable case with 𝑆
stable case with 𝑆

1.5, but is decreased in the

1.0 as in Fig. 4, where the particle flux by the resonant mode with 𝑚, 𝑛

36,14 in the inversed gradient region is shown. The unstable mode induces the particle pinch, which
increases with the larger density gradient.

Fig.3 Linear growth rate in theFig.4 Time evolutions of the particle fluxes by the ( 36,14)

reversed gradient region (𝜌 = 0.76).

mode with (a) 𝑆

1.0, and (b) 1.5.

4. Nonlinear phase results
The effects of profile relaxation and existence of multiple types of instabilities are shown in the nonlinear
phase of the calculation. Figure 5 shows the mode amplitude profiles (a) just before the nonlinear phase (t =50) and (b) during the nonlinear saturation phase (t = 120). At t = 50, the ballooning mode with 𝑛dominant, while at t = 120, a tearing mode with 𝑛

8 is

1 is dominant. The higher wavenumber components

shown in Fig. 5(a) form the ballooning structure localized on the low field side, and in the nonlinear phasethe modes spread radially to flatten the background pressure profile. This flattened profile makes theballooning modes to be stable, and the mode intensities are decreased as in Fig. 5(b). The tearing instability
grows near the rational surface with 𝑞

3, and then in the nonlinear saturation phase becomes dominant.

Next, we evaluate the particle fluxes at the radial position where the ballooning and tearing modes exist.
Figures 6(a-c) show the particle fluxes at 𝜌 = 0.6, 0.73, and 0.83, which are corresponding to the positions of
the rational surface with 𝑞

𝑚/𝑛2, 2.5 , and 3 , respectively. In Fig. 6(c), the contribution to

outward particle flux from the ballooning mode is large in the linear phase, sharply increases just in the
nonlinear saturation phase, and saturates with the relaxation of the background density profile. Figures 6(b)shows that the inward particle flux in the negative gradient region increases in the linear phase, keepincreasing after the profile relaxation in the normal gradient region. The increase even in the nonlinear phasesuggests that the contribution from the low-wavenumber components, including the tearing mode, exists inthat phase.

Fig.5 Mode structures of the excited modes (a) just before the nonlinear phase at t = 50, and (b) in the
nonlinear phase at t=120. The mode amplitudes of the density are shown.

Fig.6 Time evolutions of the particle flux at radial positions 𝜌 = (a) 0.6, (b) 0.73, (c) 0.83.

5. Conclusion

The R5F code was used to evaluate turbulent particle transport in a global simulation. The effect of thenegative density gradient was investigated by introducing a density profile with an inversed density gradientnear the plasma edge. Outward particle fluxes in the normal gradient region are dominant, and the evolutionof the density profile is largely determined by this outward transport. On the other hand, inward particle
fluxes are also observed in the inversed gradient region. The magnitude of the inward pinch is increased withthe density gradient larger than the linear threshold given by the electron drift wave mode theory. The
contribution of the inward particle fluxes was found to exist even in the nonlinear phase. The amplitude keeps
increasing in the nonlinear phase after the flattening of the pressure profile by the ballooning mode in thenormal density region.
References[1] B. Coppi and C. Spight, Phys. Rev. Lett. 41, 551 (1978).
[2] M. Yagi, et al., Proc. BPSI Meeting, 2-1 (2018).
[3] M. Yagi, et al., 26th IAEA FEC, TH/P3-21 (2016).
[4] C. T. Hsu, et al., Phys. Fluids 29, 1480 (1986).

Improvements of impurity transport modelling by integrated divertor
code SONIC and its application to predictive simulation of future
devices
Shohei YAMOTO1), Kazuo HOSHINO2), Yuki HOMMA1), Tomohide NAKANO1), Ryuichi
SANO1) and Nobuhiko HAYASHI1)
1)

QST, 2)Keio Univ.

SONIC is one of integrated edge-plasma transport codes, which consists of the fluid plasmatransport code SOLDOR, the kinetic neutral transport code NEUT2D and the kinetic impurity
transport code IMPMC. One of the most unique features of SONIC is to compute the impurity
transport processes kinetically. Implementing Multiple-Program Multiple-Data (MPMD) parallel
computing technology enables SONIC to track independently the transport of multiple impurity
species by launching IMPMC modules. Furthermore, time-dependent solver is implemented in
SONIC, so that temporal impurity transport can be analyzed. With these advanced features,
SONIC has been used for divertor operation scenario development for JT-60SA and the JA DEMO.
In this study, we have applied the time-dependent version of SONIC to the transient analysisof additional injection of Ne into Ar-seeded plasma in JT-60SA. In Ar-only case, the Ar impurities
are stagnated in top of the SOL due to the thermal force, which are the main source of impurities
into the core. It is found that only a small amount of Ne injection to the Ar seeded plasma
suppresses the Ar transport along the SOL plasma towards the top of the plasma; this is mainly
caused by the enhanced friction force due to the higher D+ parallel flow in the high-field-side
SOL towards the inner divertor. Figure 1 shows the Ar radiation power pAr and the parallel D+
flow velocity u//,D+ plotted along the field line 3.0 mm outside from the separatrix at the outer
midplane at the typical snapshots. As seen from figure 1, the acceleration of u//,D+ towards innerdivertor plasma can be seen especially in the HFS near the X point, where the Ar radiation power
was low. The initial Ne radiation takes place in the region where the Ar radiation power is low.
In summary, the possible mechanism of the screening of Ar by additional Ne injection is asfollows:(i)
The Ne radiation takes place where Ar radiation power is low at the HFS near
X-point.
(ii)
The Ne radiation cooling induces initial acceleration of u//,D+ towards innerdivertor.
(iii) The ion energy is transported towards inner divertor by accelerated u//,D+ as aconvection, which mitigates the deceleration of u//,D+ and u//,D+ become higher.
(iv) The friction force pushes Ar ions towards inner divertor by the high u//,D+. The
highly charged Ar ions also contributes to the radiation cooling at the HFS near X-point.
The Ne ions are also transported towards inner divertor and thus both Ar and Ne
impurities complements the low radiation power regions of each impurity species andthe u//,D+ is further accelerated.
(v)
More ion energy is also transported towards inner divertor by the convection
and deceleration of u//,D+ is further mitigated and leads to the decrease of Ar radiationand density at the top of SOL.
The time dependent analysis suggest that the additional Ne injection into Ar-seeded plasma is apossible channel to control impurity transport in the SOL.

Figure 1. Ar radiation power pAr and parallel D+ flow velocity u//,D+ plotted along field line 3.0mm outside from separatrix at outer midplaneSimulation of 3-D perturbation measurement in torus plasmas
using a heavy ion beam probe
R. Yoshihara1, N. Kasuya2, T. Ido2, M. Nunami3, K. Fujita3, S. Satake3
1

IGSES, Kyushu University, Kasuga, Fukuoka 816-8580, Japan

2

RIAM, Kyushu University, Kasuga, Fukuoka 816-8580, Japan
3

NIFS, Toki, Gifu 509-5292, Japan

1. Introduction
Improving plasma confinement is important to realize fusion power generation. Plasma turbulencein torus plasmas contributes to torus transport phenomena, and its three-dimensional aspect needs to
be clarified by experimental measurements. The transport phenomena are evaluated quantitatively bysynthetic diagnostics, which numerically reproduce experimental measurements using simulation data[1]. It is effective to perform numerical diagnostics on multiple types of simulation data in order toelucidate various aspects of magnetically-confined plasmas. We are developing an integrated
numerical diagnostics platform that can apply multiple analysis methods to multiple simulation data.
Here, we present numerical diagnostics of heavy ion beam probe (HIBP) for the Large Helical
Device (LHD) and the PLATO tokamak [2]. Numerical diagnostics are performed by developing HIBP
routines to three different data sets obtained by 3-D simulation code FORTEC-3D, MIPS, and R4F.

2. Heavy Ion Beam Probe
HIBP is a measurement instrument that can obtain information on the potential, density andmagnetic field inside of the plasma. Injected (primary) beam of univalent (1+) heavy ions becomes tobe divalent (2+) ions due to collisional ionization with plasma particles, and measurement of this(secondary) beam gives information at the ionization point [3]. Difference in kinetic energy between
the primary and secondary beams corresponds to the potential at the ionization point, and density
fluctuation is related to time variation of the secondary beam intensity.
3. Measurement simulation
To simulate a HIBP measurement, three-dimensional information of the magnetic field, plasma
potential, density, and temperature in the measurement space is required. 3-D field data, obtained from
FORTEC-3D [4], MIPS [5], and R4F [6] codes, are used for numerical analysis here.
The FORTEC-3D code is a drift kinetic code that generates distribution data on the magnetic flux
coordinates (𝜌，𝜃，𝜁). The MIPS code is a 3-D MHD simulation code that generates 3-D data withspatial meshes in the toroidal coordinates (𝑅，𝑍，𝜙) with equal intervals. The R4F code is a four-field
reduced MHD simulation code that uses spectral expansion in the poloidal and toroidal directions anddifferential equations in the radial direction. Thus, each simulation code handles different physicalquantities and data formats. It is necessary to use other models or assumptions, if missing physicalquantities exist. Different codes use different coordinate systems, but analysis routines to minimize
information losses are required for the synthetic diagnostic simulations.
The simulator, which simulates the HIBP measurement, calculates the trajectories of the probeheavy ions, and evaluates the potential and density including path integral effects along the trajectories.
The simulator can compare the physical quantities obtained from evaluation with and without
turbulence in 3-D magnetic field configuration. The simulator can also include the measurement error
in experiments for discussing observability of target phenomena.

4. Measurement simulation result
First, we describe simulations of the mean profile measurement, which were performed for
PLATO tokamak and LHD helical plasmas. Here, we present the results of LHD plasma by using
FORTEC-3D data. The data include a non-uniform electrostatic potential Φ on the magnetic surface
obtained from global neoclassical calculations [7]. The probe ions are monovalent gold ions Au+
injected with energy of 5.5 MeV . Figure 1 shows the ion trajectories, where the red circles represent
the ionization points. Figure 2(a) shows difference between the observed potential and the meanpotential, which corresponds to Φ . Figure 2(b) shows that the beam intensity is expected to be larger
than the instrumental noises. The beam intensity is smaller in the low field side, because the beam path
length is longer and the beam attenuation is larger in these cases.
Next, we describe the fluctuation measurement simulation. Figure 3(a) shows the time evolutionof the electrostatic potential extracted from resistive drift wave turbulence data calculated by R4F. The
calculation error is less than 1 V , and the gradual change associated with the background profilechange and the local fluctuation component with a period of about 100𝑡 are captured. Figure 3(b)shows the time evolution of the detected beam current intensity, when the ballooning instability
calculated by MIPS causes a global pressure change. In this case, the density change at the ionization
point and the time evolution of the detected beam intensity are quite different, because the density
decreases largely at the central region of the plasma and the density change along the beam path islarger. This is the case when the path integral effect cannot be ignored.

－ beam trajectory○ ionization point

Fig. 1 Heavy ion beam trajectories in the LHD configuration.

○ Calculation results△ potential data
- - mean potential

Fig. 2 Radial profiles of (a) the electrostatic potential and (b) detected beam current intensity for
HIBP simulation in LHD. The positive and negative radius in (b) corresponds to positions at low andhigh field sides compared to the magnetic axis.
－ Calculation results－ potential data

○ Calculation results○ density data

Fig. 3 Time evolutions of HIBP simulations in PLATO. (a) The electrostatic potential using R4F data
at a local point (𝜌0.37) and (b) the detected beam current intensity using MIPS data are shown.

References[1] C. Holland, et al., Phys. Plasmas 16 (2009) 052301.
[2] A. Fujisawa, AIP Conf. Proc. 1993 (2018) 020011.
[3] T. Ido, et al., Rev. Sci. Instrum. 77 (2006) 10F523.
[4] S. Satake, et al., Plasma Fusion Res. 3 (2008) S1062.
[5] Y. Todo, et al., Plasma Fusion Res. 5 (2010) S2062.
[6] M. Yagi, et al., Plasma Fusion Res. 2 (2007) 025.
[7] K. Fujita, et al., Nucl. Fusion 61 (2021) 086025.

Analysis of transport in PANTA by Multi-field SVD
T. Kodahara¹, M. Sasaki¹, Y. Kawachi², Y. Jazima¹, T. Kovayashi², T. Yamada³, H. Arakawa³
CIT, Nihon Univ.¹, NIFS², Kyusyu Univ³.
1. Introduction
Performance of magnetized plasma is dominated by turbulence driven transport. Modal decomposition
into a few degrees of freedom on real space is easier to understand detailed spatio-temporal dynamics. SVD
enables us to analyze interactions between modes.
In this study, we apply a newly proposed method to analyze the correlation among multiple physicalquantities in a data-driven manner, which is an extension of the conventional SVD method. The spatial
mode structure with a common temporal evolution can be derived by the proposed method. The method isapplied to a set of experimental data of the ion saturation current and the floating potential observed in a
linear magnetized plasma, PANTA. Turbulence driven particle transport driven is investigated with thismethod.
2. Multi-field SVD
We define the spatio-temporal matrix, which consists of the ion saturation current and the radial flow.
here, radial flow is evaluated from the floating potential with the assumption of E×B drift velocity. The
matrix is given asN(θ1 , t1 )𝑁(θ2 , t1 )⋮N(θ𝑚1 , t1 )Φij =Vr(θ1 , t1 )Vr(θ2 , t1 )⋮(Vr(θm2 , t1 )N(θ1 , t 2 )N(θ2 , t 2 )⋮N(θm1 , t 2 )

⋯
⋯
⋱
⋯

N(θ1 , t n )N(θ2 , t n )⋮N(θm1 , t n )(1)

Vr(θ1 , t 2 ) ⋯ Vr(θ1 , t n )Vr(θ2 , t 2 ) ⋯ Vr(θ2 , t n )⋮⋱⋮Vr(θm2 , t 2 ) ⋯ Vr(θm2 , t n )), where N and Vr are the ion saturation current and the radial flow, respectively. The fluctuation data ismeasured by the 64 channels probe array which is placed in the azimuthal direction at the fixed radial r=4cm
in linear plasma device, PANTA . The azimuthal position of the observation is demoted by θj (j=1～32)[1]. The SVD is applied to the matrix in order to obtain the spatial mode structures of the ion saturation
current and the radial flow with the common temporal evolution. The matrix is decomposed as
𝑚 +𝑚2

1
𝛷𝑖𝑗 = ∑𝑘=1

𝑇
𝑈𝑖𝑘 𝛴𝑘𝑘 𝑉𝑘𝑗.

(2)Here, Σ is the singular value matrix which is the diagonal matrix, and the diagonal component 𝜎1 ~𝜎𝑚1 +𝑚2
is sorted in descending order (𝜎1 ≥ 𝜎2 ≥ ⋯ 𝜎𝑚1 +𝑚2 )[2,3]. The matrix U corresponds to the spatial structure.
The matrix ΣVT has the temporal evolution of turbulence common to the ion saturation current and theradial flow. the above method is based on orthogonal bases in a data-driven manner.
3. Transport analysis by Multi-field SVD

Based on the Eq. (1), fluctuation data of the ion saturation current and the flow are combined in spatial
direction, and the singular value decomposition (SVD) is applied to the combined observation matrix. Thisis the Multi-field SVD, where the common characteristics can be deduced between the different physicalquantities. The spatial structures Uik with common temporal evolution ΣkkVTkj of fluctuation of the density
and the radial flow are evaluated, are shown Fig. 1. The phase of both spatial structures can be systematicallyevaluated. This phase difference is appeared in disposition of typical resistive drift wave [4].
Mode α of the density and mode β of the radial flow obtained from results of Multi-field SVD isrespectively defined as Nα and Vrβ, and transport is evaluated by the following̅̅̅̅̅̅̅̅< 𝛤̅𝑟 >= ∑𝛼,𝛽 < 𝑁
𝛼 𝑉𝑟𝛽 >

(3)

In addition, the results of conventional SVD is applied to Eq. (3) to evaluate transport. In the Multi-field
SVD and the conventional SVD, the transport which the combination of the such mode drives is shown in
the Fig. 2. Transport using results of Multi-field SVD is evaluated at only combinations of mode satisfying
α = β . Transport using result of conventional SVD is complicated because transport is caused by
combinations of various mode. It can be said important to simultaneously analyze fluctuation of the density
and the radial flow because of stemming from the mode orthogonality of the SVD. Multi-field SVD alsoenable to greatly reduce the number of combinations of mode.

Fig. 1 (Left panel) The spatial structure of

Fig. 2 The transport which the combination

mode1,3,5. (Right panel) The common temporalof the such mode drives, using results

evolution of N and Vr of mode 1,3,5.

obtained from the Multi-field SVD (leftpanel) and the conventional SVD (rightpanel)

References[1] T. Yamada, et al, Nature Physics, 4.9, 721-725, (2008)[2] M. Sasaki, et al, Plasma Physics and Controlled Fusion, 63.2, 025004 (2020)[3] Y. Asahi, et al, Physics of Plasma, 28.1, 012304, (2021)[4] W. Horton, et al, Reviews of Modern physics, 26.4, 042306, (2019)