J. Plasma Fusion Res. Vol.100, No.7 (2024) 301-303

小特集

古くて新しい新古典輸送理論の新展開
Recent Progress of Old yet New Neoclassical Transport Theory

1 ．はじめに
1. Introduction
佐竹真介

SATAKE Shinsuke
核融合科学研究所

構造形成・持続性ユニット（原稿受付：2024年 2 月21日）

磁場閉じ込め核融合の研究は半世紀以上に渡って続け

そのような核融合理論研究の歴史の中で，新古典輸送

られており，最近では JT-60SA のファーストプラズマ成理論はトーラスプラズマ中の荷電粒子の衝突拡散現象を功やフランスにおける ITER 装置の建設など，大型超電扱う理論として比較的早期から研究が進められてきてお

導コイル装置における DT プラズマの定常燃焼の実現にり，1980 年頃には現代でもよく引用される新古典輸送理向け，一歩ずつではあるが着実に研究が進められている．

論の体系がまとめられた［1-4］．そのような「ベーシック」

またこの数年間で見ると，海外だけでなく日本においてな新古典輸送理論は現在まで続く核融合研究の基盤とも民間資本による多くのスタートアップ企業が核融合炉なっており，その中では厳密な方程式系の導出も示され実現に向けた研究に参入してきており，トカマク型，ヘていたが，それを実際に応用する際には，計算機が未発

リカル型の両方で新しいコンセプトの装置の建設とそれ達であったため定量的な正確性よりむしろ定性的な理解による核融合発電の実証に向けた実験がこの先多く進めに重点が置かれていたと言っても過言ではないであろう．

られようとしている．しかしこの現在の活況は数年前は新古典輸送計算を実際に行うために採用された多くのモ予想もつかなかったことであり，核融合研究の歴史におデルは，計算する上での簡便性のために様々な近似や仮いて，新しい実験装置を建設することは多くのコストと定のもとに成り立っている．当時の計算機の性能を考え時間を要することから実行に移すハードルが高く，これればこれは至極当然のことであり，筆者はこれを批判すが核融合研究に長い年月を要している原因の一つと筆者るわけではないが，そうした簡易モデルに基づく理論がは考える．そのような状況においても，核融合プラズマ

標準的な新古典輸送理論として教科書等で紹介され，まに関する理論研究においては 21 世紀に入って以降各段のたそのハンドリングのしやすさから実験解析においても進展があったと言え，その大きな要因は計算機の発達に多く採用される状況が長年続いたことによって，核融合伴う数値シミュレーション法の高度化にある．今では核研究者の中に「新古典輸送理論は既に完結しており，し融合プラズマの理論研究において単に従来より大規模・かも容易に計算可能」という見方が，特に新古典輸送が高解像度なシミュレーションを行うだけではなく，それ低く抑えられ，かつ磁場配位が単純な故に簡易的な計算までの計算法で無視されてきた効果を取り入れたより現法でも精度のよい新古典輸送計算ができるトカマクの研実に近いモデルに基づくシミュレーション法を開発した

究コミュニティを中心に広まってしまっているように思り，個別に扱われていた現象を統一的に扱うマルチスケー

われる．

ル，マルチフィジックス（multi-scale, multi-physics）モ

これと比較して，2000 年以降の理論研究で特に目立っ

デルのシミュレーションを行ったりすることで，これまた進展をしているのがジャイロ運動論に基づく微視的乱でに見過ごされていた現象，これまでの計算法では再現

流輸送の研究である．ジャイロ運動論による乱流輸送のできなかった現象を数値シミュレーションで研究するこ研究も，その初期には平板（スラブ）配位を使いプラズとは当然のように行われるようになってきている．

マのごく一部を切り出した簡易モデルによる研究から始National Institute for Fusion Science, Toki, GIFU 509-5292, Japanauthor’s e-mail: satake.shinsuke@nifs.ac.jp

301Ⓒ2024 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research

まった．しかしその後，スーパーコンピュータの発達とVol.100, No.7

July 2024

外からは意識されていないこともあるのが現状である．

共にジャイロ運動論の研究は進展を続け，実際的な 3 次

筆者は新古典輸送の理論シミュレーション研究が 2000

元トーラスプラズマの磁場配位を入れ，イオンと電子，

年以降，教科書的なよく知られているものから，計算機

さらに不純物イオンを含め，粒子種毎に異なるスケールとシミュレーション技法の発達と共により高度なものへの乱流の相互作用までを同時に解くレベルにまでシミュと進展してきた経緯を自らも携わりながら追ってきた．

レーションの解像度と規模が進展している．また，そのそのような立場から，ここで述べてきた新古典輸送理論・

進展の過程においてジャイロ運動論のより数学的・物理

シミュレーションに対するある種不幸な誤解を解き，新的厳密性を備えた定式化と，そのモデルに基づいた大規

古典輸送理論そのものを研究している研究者以外の大多模数値計算を効率的かつ高精度に行うための様々な数値

数の利用者としての立場の核融合研究者が，予測研究や計算技法の開発も絶え間なく進められてきた．このよう実験解析等に新古典輸送計算を利用する際に，その用途な経緯から，核融合コミュニティの中でジャイロ運動論に適したモデルの計算法を選択し利用できるようになっによるプラズマ乱流研究を完成されたものとみなす研究て欲しいと考えてこの小特集を執筆することにした．ま者はまずいないであろう．実際，シミュレーション研究

た，新古典輸送理論は今でもまだ進化し続けており，最の進展とともにトーラスプラズマ中の微視的乱流輸送現

近どのような拡張の取り組みがなされてきているかにつ象に関する新しい知見の数々が産まれ続けており，さらいてもいくつか紹介したい．こうした新しい拡張は，定にこれからはそうした知見を新装置の設計や運転予測・

常核燃焼プラズマ中の輸送を評価する上で重要になると制御に活かそうという取り組みも始まっている．

思われる物理を新古典輸送計算に取り入れるために必要

それでは，そのようにジャイロ運動論シミュレーショ

であり，この小特集をきっかけに我こそはと思う若手にンの研究が進展している間，新古典輸送現象の理論シミュ

ぜひ更に発展させてもらいたいと願っている．なお，こレーション研究は止まってしまっていたのかと言うと決こに述べたシミュレーション法に対する古い観念に基づしてそのような事はない．特にヘリカル系プラズマにおく誤解や，計算に採用されているモデル固有の特徴や近いては，新古典輸送現象は乱流輸送同様にプラズマ閉じ

似による限界などに対する理解不足の問題は新古典輸送

込め性能を評価する上での最重要指標の一つであるため，に限った話ではない．MHD 平衡・安定性コード，乱流シ

2000 年以降も従来の新古典輸送計算モデルで使われていミュレーションコードなど，現在の核融合研究においてた近似を排したより定量的精度の高い計算法の開発が，ヘ

複数の，他人が開発したコードを組み合わせて利用する

リカルプラズマの研究分野を中心に 2000 年以降も進めら

ことは理論研究者にとっては当然のことになっているが，れてきた．ここでも，コンピュータの性能と数値計算技法，

そのそれぞれの計算法の中でどのような近似が使われ，両方の発達とともに以前より現実的な仮定に基づく計算

どの程度の定量的，定性的な正確さが期待できるかをしっ

モデルが新古典輸送計算に適用されている．さらに，従かりと理解して使っているかどうか，この小特集をきっ

来は微視的乱流を扱うジャイロ運動論と，新古典輸送現かけに一度考え直す機会になれば幸いである．

象を扱うドリフト運動論はスケール分離可能という仮定

本小特集の構成は以下の通りである．まず第 2 章では，のもとでそれぞれ別個に計算されてきたが，この 10 年間

新古典輸送とはそもそも何なのか，電磁流体力学（MHD）

でジャイロ運動論シミュレーションの中に新古典輸送の乱流や古典輸送，微視的乱流との区別について，特にジャ効果を含めた計算ができるようになってきた．乱流現象

イロ運動論とドリフト運動論の導出過程を概観しながらは基本的に非線形なモード間のカップリングが重要であ定義していく．続いて新古典輸送現象について，特に重り，それゆえに実際に乱流輸送のレベルを評価するには要となる小半径方向の粒子・熱フラックスやブートスト

数値シミュレーションに頼らざるを得ないのと比較して，ラップ電流と，新古典粘性，Coulomb 衝突による粒子種

新古典輸送現象は磁場配位とプラズマの分布を与えれば間摩擦などがどのような関係式で結ばれているかを見る．

基本的に確定的に解が求まる類の問題である．このこと

また，トカマクや準軸対称系と，大型ヘリカル装置（LHD）

も乱流現象と比較して新古典輸送の計算は容易であるとや Heliotron-J,W7-X などの非軸対称系配位との大きな違時に勘違いされる要因になっているかも知れないが，現いである両極性径電場についても説明する．なお，本小在の新古典輸送計算は採用するモデルの近似の度合いに特集ではよく教科書レベルの新古典輸送理論で紹介され応じて，手計算レベルのものから，デスクトップパソコ

る，トカマクのバナナ - プラトー - Pfirsch-Schlüter 領域

ン 1 台程度のもの，さらにジャイロ運動論シミュレーショ

など衝突領域ごとの新古典輸送が具体的にどのようになンと同様にスーパーコンピュータを使った大規模並列計

るか，ヘリカル系ではどうか，などの定性的な説明は行算が必要なものまで，様々なレベルのものが存在していわない．ここでは，新古典輸送現象を扱う方程式系がどる．計算モデルの選択によってはある種の現象を評価でのような定式化によって導かれるかについて理解してもきない，あるいは過大・過少評価してしまうことがあり，らうことを主な目的とする．次に第 3 章では，ドリフト理論予測や実験解析に新古典輸送シミュレーションを用運動論を実際にシミュレーションで解く際に，どのよういる際にはそれが採用している計算モデルの特徴をよくな近似が用いられ，それぞれの近似がどのような影響を理解した上で利用する必要がある．しかし，このような計算結果に与えるかについて，いくつかの主要な新古典

詳細についてはコードの開発者や新古典輸送の専門家以

輸送コードを挙げながら説明する．さらに，新古典輸送
302

Special Topic Article

1. Introduction	

S. Satake

計算の更なる拡張の方向性として，いくつかの研究を紹

この小特集をきっかけにより深く理解していただけると介する．これらは，中性粒子ビーム入射（Neutral beam

幸いである．また，研究の中で既に新古典輸送コードを

injection, NBI）などの外部加熱や核融合反応による α 粒

ユーザーとして利用している方，あるいはこれから利用

子による加熱の効果など，非熱平衡的な分布関数に対すする方は，3．1 章や 5 章の議論から新古典輸送コードが採る取り扱いと，多粒子種プラズマの新古典輸送シミュレー

用する計算モデルに含まれる近似や制約の概観がつかめ，

ションによる不純物イオンの新古典輸送現象に関する研さらに各章を読むことで理解が深められると思うので，究である．さらに 4 章では，新古典輸送シミュレーショ

本小特集が新古典輸送現象に関する研究を進める上でのン研究の最先端として，ジャイロ運動論による微視的乱一助となることを願う．

流と新古典輸送現象の相互作用に関する研究例を紹介してもらう．最後に 5 章で本小特集の内容に関するまとめ参考文献を述べる． 4 章で紹介するような最新のシミュレーション［ 1 ］F.L. Hinton and R.D. Hazeltine, Rev. Mod. Phys. 2,239 (1976).
［ 2 ］S.P. Hirshman and D.J. Sigmar, Nucl. Fusion 21,1079 (1981).
［ 3 ］L.M. Kovrizhnyikh, Nucl. Fusion 24, 851 (1984).
［ 4 ］K.C. Shaing and J.D. Callen, Phys. Fluids 26, 3315
(1983).

モデルを用いた研究においても，その中で新古典輸送（ドリフト運動論）がどのようなモデルに基づいて数値的に解かれているかを理解することは，計算モデルの妥当性を説明する上でも非常に重要な点であるので，ジャイロ運動論を主に研究している研究者にも，カウンターパー
トとしての新古典輸送計算がどのように行われているか，

303

J. Plasma Fusion Res. Vol.100, No.7 (2024) 304-313

小特集 古くて新しい新古典輸送理論の新展開

2 ．新古典輸送理論の概要
2. Overview of Neoclassical Transport Theory
佐竹真介

SATAKE Shinsuke
核融合科学研究所

構造形成・持続性ユニット（原稿受付：2024年 2 月21日）

新古典輸送理論の礎となるドリフト運動論方程式の導出から，新古典輸送方程式の導出と，モーメンタムバ
ランスや両極性径電場の決定など，新古典輸送現象において重要となる事柄についてまとめる．また，ジャイロ
運動論との違いについても概観する．
Keywords:
Neoclassical transport theory, drift-kinetic equation, gyrokinetic theory, numerical simulation

2．1

輸送現象のクラス分け

われるのは閉じ込め磁場と電場（B，E）の他には密度 n，温トーラス型のプラズマ閉じ込め装置は，トーラスを周

度 T，圧力テンソル P，流体速度 V といった一般の流体力

回する磁力線によって作る入れ子状の磁気面構造によっ学でも使われる量であり，後に述べる運動論的モデルでて，プラズマを真空容器壁に接触しないように隔離し，扱う速度空間分布の情報を必要としない（Alfvén 共鳴モー高温のプラズマを長時間安定に保持することを目的とす

ドを駆動する高速粒子の分布など，速度分布の情報が重るものである．プラズマを構成する電子とイオンは磁力要な現象も一部含まれる）．また，一様磁場中でジャイロ

線に巻きつくジャイロ運動をしながら磁力線方向に走り運動する荷電粒子同士の Coulomb 衝突による，Larmor 半

回る．従って端のないトーラス状の磁気面を構成するこ径をステップサイズとする拡散現象を古典拡散と呼ぶが，とで荷電粒子を長時間閉じ込める，というのが磁場閉じ後で見るようにこれは流体力学的変数によって記述可能

込め核融合炉の基本的なアイデアである．しかしながらなため，流体的な輸送現象の 1 つと見なすこともできる．

荷電粒子は永久にトーラス中を周回し続けるわけではな

MHD 不安定性に伴う輸送は核融合炉の定常運転を実現す

く，様々な原因によって磁気面を横切る方向に拡散される上で十分に抑制されるように磁場配位や運転条件を設る．また正負の荷電粒子の流速の差はプラズマ中を流れ計・調整しなければならない．また，古典拡散による熱る電流となり，プラズマ自身を閉じこめている磁場を変輸送のタイムスケールは核融合炉に求められる熱閉じ込化させる．このような閉じ込め磁場中の粒子や熱の流れ

め時間に比べて十分に遅いため，流体的な輸送現象だけを扱うのが磁場閉じ込め核融合の研究分野における輸送を想定するならば，MHD 安定なプラズマは核融合炉の炉理論（transport theory）である．

心条件を容易に満たす．しかしながら実際には，古典拡
散よりずっと大きい輸送をもたらす運動論的な輸送現象

磁場閉じ込めプラズマ中で起こる輸送現象には，大別すによってプラズマの閉じ込め性能は決定されている．

ると流体的なものと運動論的なものの 2 つがある．前者
は， 電 磁 流 体 力 学（Magneto-Hydro-Dynamics，MHD）

運動論的な輸送現象は，安定的な MHD 平衡磁場が存理論で扱われる閉じ込め磁場の巨視的な流体力学的不安在しているという前提のもとで，荷電粒子の平均自由行

定性に伴い，閉じ込め磁場そのものがマクロ（スケール

程 lmfp が閉じ込め磁場のスケール長（～大半径 R）程度か

長 L ～小半径 a）かつ高速（時間スケール t ～ L/vth，i，vth，i

それ以上となる，高温で低衝突性な磁場閉じ込めプラズ

はイオン熱速度）に変化し，時にプラズマの圧力分布自マにおいて特徴的に現れる現象である．平均自由行程が体が大きく変化するような巨視的な粒子・熱の吐き出し

長く衝突による緩和が効きにくいこと，ジャイロ運動を現象である．この代表例としては，キンクモード，バルーする粒子の旋回中心（ガイディングセンター）がトーラ

ニングモード，Alfvén 共鳴モードなどの MHD 不安定性

ス磁場の不均一性や電場の影響で磁力線を横切る方向にが挙げられる．このようなマクロな不安定性に伴う輸送

ドリフト運動すること，そのドリフト軌道が粒子の磁力では，磁場に閉じ込められたプラズマのエネルギーを大線平行及び垂直方向の比（v‖ : v⊥）に応じて様々に変化す

きく失うことになる．MHD 理論で現象を記述する際に使ること，プラズマの温度密度が磁気面関数になっており

National Institute for Fusion Science, Toki, GIFU 509-5292, Japanauthor’s e-mail: satake.shinsuke@nifs.ac.jp

304Ⓒ2024 The Japan Society of Plasma
Science and Nuclear Fusion Research

Special Topic Article

2. Overview of Neoclassical Transport Theory	S. Satake

磁気面を横切る方向に温度密度勾配が存在すること，こつ［5］．このように，新古典輸送は磁場閉じ込め核融合装

れらのことの帰結として速度空間における荷電粒子の分置の設計や閉じ込め性能予測にとって重要な要素の一つ

布関数には熱平衡状態である局所 Maxwell 分布からのずであり，核融合研究の最初期から今日に至るまで継続的れ（非等方性）が残る．この速度空間の非等方性から，に研究が続けられている．

一つには圧力テンソル P の磁力線方向成分と垂直成分に2．2

関する非等方性と空間的非一様性を生み，その応力（プ リフト運動論方程式とジャイロ運動論方
ド
程式

ラズマの運動論では一般に粘性と呼ばれる）粘性 ∇・P と
Coulomb 衝突による摩擦力のバランス関係が運動論的な

ここからはまず，磁場閉じ込めプラズマの運動論的記

輸送現象を引き起こす．これを新古典輸送現象と呼ぶ．

述から，新古典輸送現象を扱うドリフト運動論方程式と，何をもって「新」古典なのかと言えば，今説明した現象

微視的乱流輸送を扱うジャイロ運動論方程式がどのようはトーラス磁場が非一様で曲がっているために生じるドに導出され，どのように両者が分離されるかの概要を見

リフト運動に起因するためである．これに対して古典拡ていく．歴史的に両者はスケール分離が可能であるとい散現象は一様な磁場中でも起こる，ジャイロ運動をする

う前提のもと，それぞれ個別の現象として扱われてきた

荷電粒子の衝突拡散現象を指す．運動論的な輸送現象にが，最近ではスーパーコンピュータの発達に伴い両者をは新古典輸送の他に，Larmor 半径程度の微視的な電磁場

分離せずに 1 つの運動論的モデルとして数値計算するシの揺らぎによって生じる波（ドリフト波）によって生じミュレーションも発展してきている．そこではこれまで

る乱流現象がある．これは位相空間においてある種の共想定されていなかった乱流と新古典輸送の相関関係があ鳴を引き起こし，磁気面を横切る方向に微視的な渦を成ることが発見され始めている．本節ではまず，従来から長させ，プラズマの粒子・熱閉じ込め性能を劣化させる．

の定式化に従って 2 つの運動論方程式の導出過程を概観

これを異常輸送，あるいは（微視的）乱流輸送などと呼ぶ．

する．

新古典輸送，乱流輸送ともに磁気面を横切る方向の温度

運動論的記述のベースとなるのは，空間 3 次元 X，速密度勾配を駆動力として生じる．運動論的な輸送現象は度空間 3 次元（U，μ，ξ）の 6 次元位相空間に荷電粒子のおMHD 的なものに比べると穏やかかつゆっくりした時間ス

けるジャイロ中心の分布関数 F（X，U，μ，ξ，t）である．ここケールで起こるが，核融合炉の定常的な閉じ込め性能をで，X はジャイロセンター位置，U はジャイロセンターの評価する上で重要な役割を果たし，実験においても実質磁力線方向速度，μ は磁気モーメント，ξ は磁力線周りの的にプラズマの熱閉じ込め時間を決定している．

ジャイロ運動の位相角，t は時間を表す．以下，（X，U，μ，ξ）

新古典輸送と乱流輸送の大きさを比べると，新古典輸をまとめて Z と表す．これらの変数の厳密な定義の説明

送は一般的にドリフト軌道の磁気面を横切る幅に応じては文献［6，7］に譲るが，近似的には荷電粒子の位置と速度大きくなるが，このドリフト軌道幅は磁場配位にある種

z＝（x，v）を 用 い て X＝x－ρ ＝ x－mB × v⊥/eB2，U＝v‖，の空間的な対称性を持たせることによって小さくできる．

2/2B である．なお，文献［6，7］においてはジャイ
μ＝mv⊥

そのため核融合実験装置では磁場形状を工夫することで，ロ平均された粒子位置のラグランジアンが ξ に依らないよ

程度の差こそあれ新古典輸送は低いレベルに抑えられてうに定義された「ジャイロ中心」と，初学者向けの教科

いる．軸対称磁場を持つトカマク装置はその最たるもの

書で一般的に使われる，粒子位置から Larmor 半径 ρ だけであるが，他にも HSX，W7-X，CFQS などの先進的なス

移動したガイディングセンターとを厳密に使い分けていテラレータ装置は，複雑な形状のコイルによって新古典るが，本章では細かい区別には立ち入らずジャイロ中心

輸送を抑制する磁場配位を作ることが主要な設計指針とと呼ぶことにする．今，新古典輸送や乱流輸送で扱う輸なっている．一方，微視的乱流はプラズマの圧力勾配が送現象の時間スケールがジャイロ運動の時間スケール～

強くなるほど一般的に不安定になり，高温プラズマでは

Ω－1 より十分長いと仮定すると（Ω ＝ eB/m はジャイロ周

微視的乱流によって閉じ込め性能が決まっていることが波数），μ は断熱不変量と見なせる．また，背景の電磁場

これまでの様々な装置での実験から評価した閉じ込め時（B0，E0）やプラズマの温度密度（n0，T0）も同様にゆっくり間のスケーリング則から示唆されている［1，
2］
．しかし，

時間変化しており，それらの空間の勾配スケール L に比べ

例えば新古典輸送が比較的大きい LHD では新古典輸送とLarmor 半径 ρ は十分短いと仮定する．ドリフトオーダリ

乱流輸送のレベルが同程度になる場合も見られ［3，4］，まングと呼ばれるこの空間スケールに対する分離を，パラた 4 章で紹介するように乱流と新古典の輸送の釣り合い

メータ δ ＝ ρ/L≪1 を用いて表現する．

がプラズマの不純物輸送レベルを決定しているとみられジャイロ中心分布関数 F の時間発展は以下の運動論方

るケースも見つかっている．軸対称なトカマク装置では程式で記述される（以下の導出過程については［8］を参新古典輸送による磁気面を横切る径方向の輸送は乱流に照）．

比べ十分小さいが，ブートストラップ電流と呼ばれる，  ∂ ⎫⎧ ∂ ̇ ∂⎨ ＋X⎬ 〉ξ＝〈C 〉ξ
＋U̇
〈F
⎩ ∂t
∂X
∂U ⎭

径方向の圧力勾配によって新古典的なメカニズムで駆動される磁力線方向の自発的な電流は，閉じ込め磁場の形

ここで Ȧ ＝成にプラズマ内部のトロイダル電流を必要とするトカマ
クにおいて外部駆動電流の割合を下げる重要な役割を持
305

（2．1）

dA
〈⋯〉
，
ξ はジャイロ角 ξ に対する平均を表し，dt

Journal of Plasma and Fusion Research

右辺は Coulomb 衝突項である．また μ̇ ＝ 0 を利用している．ジャイロ中心の運動方程式 Ż は粒子運動のジャイロ

Vol.100, No.7

July 2024（  2 ＋ μB ）となること，またこれの磁気面平均を取り
maU2

平均から求められ（厳密には摂動場の効果を含めた正準0変換を用いてジャイロ中心及びその運動論方程式は導か

準中性条件を課すことで背景密度温度（n0a，T0a）及びアンれる），以下のように与えられる（本文中では原論文［6-8］

サンブル平均した静電ポテンシャル Φ0 ＝〈ϕ〉ens が磁気面と異なり MKSA 単位系で表記する）．

（

関数となることが示される．なお，上式とこれ以降は，） （

必要に応じて粒子種の添え字 a，b を記す．次に（2．1）式の）e ∂〈ψ〉
∂A*
1
μ
ξ
̇
X＝
U＋    B *＋b× ∇B0＋∇
〈ψ〉
，
ξ＋
B＊
e
m
∂U
∂t
‖

＊

（

δ1 次の式を考える．
（2．1）式のアンサンブル平均を取ると，
fˉa1 ＝〈Fa1〉ens に対する次の運動論方程式を得る．

（2．2a）

）

B
∂A*
̇
U＝－  
〈ψ〉

・μ∇B0＋e ∇ξ＋   ，
mB＊∂t‖（2．2b）e2 ∂〈ψ〉ξ̇
ξ＝Ω＋     ．m ∂μ（2．2c）

Ub・∇－

（

）

eaE‖U
L ˉ
＋ vda・∇－    Fa0＝∑ Cab（fa1）．
T0a
b

（2．4）

L
は線形化衝突オペレータで，
ここで右辺の Cab

なお上式において，背景磁場を平衡部と乱流による揺動 のように分け，これと静部分，B ＝ B0＋B＝∇×（A0＋A）L ˉ
Cab（fa1）＝Cab（fˉa1，Fb0）＋Cab（Fa0，fˉb1）

電ポテンシャル ϕ（x，t）から以下のジャイロ中心座標の関数を導入した．なお，b＝B0/B0 であり，（2．2）式は ξ に依（2．5）

を略記したものであり，衝突相手の粒子種の和 ∑b には自己衝突 b ＝ a も含まれることに注意．また，vda は以下の存しないことに注意．

ような磁力線を横切るドリフト速度を表す．

A（X＋ρ，t），
ψ（Z，t）＝ϕ（X＋ρ，t）－v（Z，t）
・
A＊
t）＋
（Z，t）＝A（X，
0

μ
∂ ˉ
fa1
b・∇B0
ma
∂U

1
vda＝  b×（maU 2b・∇b＋μ∇B0＋ea∇Φ0） （2．6）
eaB0

m
Ub，B＊
（Z，t）＝∇×A＊，e＊
B＊‖ ＝B ・∇b．

（2．3）

これは左から順に，曲率ドリフト，grad-B ドリフト，お次に，ジャイロ中心の分布関数 F（Z，t）を背景電磁場の揺動に対する応答として現れる部分 Fˆ と，統計平均（アよび E×B ドリフトと呼ばれるもので，無摂動磁場と静ˆンサンブル平均）を取ったもの〈F 〉
ens に分ける．F の特揺動と無関係である．また，E‖＝－b・［∇Φ1 ＋ ∂A0∂t］は電ポテンシャル B0 および Φ0 から与えられており，場の徴的な周波数は ω ～ vt /L であり，前述の仮定より ω≪Ω

磁力線方向の電場を表しており，磁気面上のポテンシャ

り，その特徴的な周波数は δ ω と仮定される（いわゆる輸ル非一様性は Φ1～ δΦ0 のオーダーと仮定される．（2．4）式
をドリフト運動論方程式と呼び，その解 fˉa1 が生み出すである．一方，〈F 〉ens は Fˆ より十分ゆっくり変化してお2送（transport）オーダリング）．Fˆ の空間スケールは非等方性を持ち，磁力線垂直方向には Larmor 半径程度の微視

トーラスプラズマ中の輸送を新古典輸送と呼ぶ．一方，Fˆ a1 ＝〈Fa1〉ens－fˉa1 に対する運動論方程式は（2．1）式の δ1 次

的擾乱を含んでおりそのスケール長は ρ 程度，磁力線平行の式から（2．4）式を引くことで導かれる．

方向には L 程度と仮定される．また，揺動成分の振幅に
ˆ〈F 〉ens ∼ δ を仮定する．なお補足すると，アンサ
関して F/

∂
μ
∂ ˆ
Fa1
da）
＋（Ub＋vda＋v・∇－ b・∇B0
ma∂t∂U

ンブル平均とは統計力学における概念で，対象としてい

eaFa0
μ
∂
b・∇B0
＋   （Ub＋vda）
・∇－〈ψˆ a〉ξ
ma
T0a∂U

る系と巨視的には同等な多数の仮想的な系における，様々な微視的状態に対しその出現確率の重みをかけて平均をL
da・∇Fa0＝∑ 〈Cab＋v
（Fa1）〉ξ

取ったものである．エルゴード仮説に基づけば，ジャイ
ロ運動論において上述の揺動成分 Fˆ の波長より長くかつ

b

（2．7）

流揺動の時間スケールより長くかつ平衡量の時間変化よˆt）＝ϕ（X＋ρt）－da ＝ b × ∇〈ψˆ a〉ξ /B0（ψˆ（Z，こ こ で，v
aa，v（Z，t）
A（X＋ρa，t））は摂動電磁場によるジャイロ中心・

り短い時間で取った時間平均はアンサンブル平均と同等のドリフト速度成分であり，乱流場による軌道の乱れをであると解釈される．F のドリフトオーダリング δ による表している．（2．7）式をジャイロ運動論方程式と呼び，その解 Fˆ a1 が生み出す輸送が微視的乱流や異常輸送などと呼平衡量の勾配長より短いスケールの局所空間平均や，乱展開，F ＝ F0 ＋ δF1 ＋…を考え（2．1）式に入れると，最低次（δ－1 次）及びその次のオーダーの式から F0 ＝〈F0〉ens 及L
ばれるものになる．なお，線形化 Coulomb 衝突項 Cab
はび F1 はジャイロ角に依存しないことが示される．特に δ0fv）
本来，粒子位置と速度（x，v）で表現された分布関数（x，

次のオーダーの式に対するオペレータである．（2．4）式，（2．7）式において，
L
Cabの引数として現れる fˉa1，Fˆ a1 は共にジャイロ角 ξ に非1∂
Ub・∇－ b・
（Fa0，Fb0）（ μ∇B0＋ea∇Φ0） Fa0＝∑ Cabma∂Ub
か ら，Fa0 が Maxwell 分 布 Fa0 ＝ na0

（

t）を Z
依存であるが，変数変換 z ＝ z（Z，t）を用いて Fˆ（z，1で表わしたものと解釈される．そのため，（2．7）
式では衝）ma 3/2
1
exp －
  
T0a
2πT0a

突項自体のジャイロ平均を取っている．一方，（2．4）式の306

Special Topic Article

2. Overview of Neoclassical Transport Theory	S. Satakeて，V ′ ＝ dV/ds ＝ dθdζ gB（ gB は磁気座標系のヤコドリフト運動論方程式では，衝突項のジャイロ角依存性に関係する部分（有限 Larmor 半径効果）は 2．3 節で扱う

ビアン）である．

古典輸送の導出で分離され，（2．4）式以降の定式化におい

一方，上記のような粒子フラックスの分類は，以下のては無視される．

fv）を用いたモーメンタム
ように粒子位置の分布関数（x，バランスの式からも導くことができ，こちらの方が新古

2．3

 ラックスの定義から見た新古典輸送と乱フ
流輸送

典輸送を生み出す機構を理解しやすい．まず，前節と同様に電磁場をアンサンブル平均成分（B0，E0）と揺動成分
ˆ E）
ˆ に分ける．同様に，分布関数も f ＝ f0 ＋ fˆと分離する．
（B，

入れ子状の閉じ込め磁気面が存在しているトーラスプ
ラズマで重要な輸送過程は，磁気面を横切る粒子及び熱

粒子分布関数の運動論方程式

フラックスである．今，磁気面ラベル，ポロイダル角，

（

）

df
e⎧
∂
1
⎫ ∂
f＝C（ f ）＝＋v・∇＋ ⎨E＋ v×B⎬・
⎭ ∂v
c
dt
m⎩
∂t

トロイダル角（s，θ，ζ）を用いて磁気座標系（Boozer 座標系［9］）が構築されており，平衡磁場がにこの分離を適用し，アンサンブル平均を取ると

B0＝∇s×∇θ＋ι∇ζ×∇s＝G∇ζ＋I∇θ＋β＊∇s（ ∂t＋v・∇＋ m {E ＋v×B }・∂v ）f ＝〈C（f）〉 ＋D，e∂で与えられるものとする．ただし s は磁気面に囲まれるト
ロイダル磁束を 2π で割ったもの，ι（s）は磁力線の回転変I s）は磁気面 s に囲まれたポロイダ
換を表し，G（s）および （

0

0

〈

∂

0

（2．11）

〉

ˆ
e 
・∂f
D≡－
{E＋v×B}m∂v ens

ル及びトロイダル電流束を 2π で割ったものである．β＊はPfirsch-Schlüter 電流に関係する部分であるが，低ベータ

ens（2．12）

プラズマでは無視できる．新古典フラックスの表式を得を得る．D は乱流場と分布関数の揺動成分の相関を表してるために，最初のジャイロ運動論方程式（2．1）式に立ち返いる．次に，（2．11）× mv として速度空間でモーメントをり，これの速度モーメントを取ることを考える．ジャイ

取ると，次のモーメンタムバランスの式となる［11，12］．

ロ中心座標の速度変数（U，μ，ξ）に対して，次の積分を a 種∂（mnu）＝－∇・P＋en（E0＋u×B0）＋R1＋K1．（2．13）∂t

粒子種の速度空間積分と定義する．
B＊‖ad3vc≡ dξ dU d μ  ．
ma



  

こ こ で，n ＝ d3vf0，nu ＝ d3vvf0，P ＝ md3vvvf0，（2．8）

R1 ＝ md3vv〈C（ f ）〉ens，K1 ＝ md3vvD である．さらにこ

ここで，B＊‖a/ma は z から Z への座標変換のヤコビアンであの式に B0×∇s を内積して磁気面平均を取ると，径方向フる．
（2．1）式を（2．
8）式で積分しアンサンブル平均をとる

ラックスの表式を得る．

と，次式の粒子保存の式（連続の式）を得る［10］．
∂nac
  ＋∇・
（ΓaNC＋Γa＋Γac）＝0，
∂t

d v F ，Γ ＝d v （Ub＋v ）fˉ ，Γ ＝d v 〈vˆ Fˆ 〉 ．

nac＝

3

NC
a

a

c a

3

3

c

c

da

da a1 ens

a1

（2．9b）〈〉〈〉
R ×B ・∇s
K ×B ・∇s∂ n      ＋
b×u ・∇s ．
＋〈     〉
＋
〈
〉
〈
） 〉
eB
eB
∂t（ Ω


（2．9c）

ここで，2．2 節の結果からドリフトオーダリングの δ1 まで

B0×∇s・∇・P
nE0×B0・∇sΓs≡〈nu・∇s〉＝－       ＋      eB2B2

（2．9a）

1

0
2

1

0
2

（2．14）

で粒子分布関数 f が次のように展開できる［8］．
（2．9d）
ˆ （X＝x－ρ，
fv）＝f0＋f＝Fv）＋fˉ1＋fˆ（x，0nagc はジャイロ中心密度，ΓaNC および Γa はそれぞれ新古典と乱流の粒子フラックスを表す．Γac は古典拡散に対応し，

ˆ
eϕ（x，
t）
ˆ
v）－F0 ρ・S0＋fˉ1－    F0＋h，
＝F（x，
0
T0

これは（2．1）式右辺の衝突項の中に入る分布関数が粒子位

（

）

mv2 5
∇p0 e
S0≡  ＋∇Φ0＋   － ∇ln T0，
p0
T0
2T0 2

f
置の分布関数（z）
ではなくジャイロ中心分布関数 F（Z）であることからジャイロ平均を取った時に残る部分があり，

ˆe〈ψ（X，v，t）〉ξ
ˆv，t）＋       F0h（X，v，t）＝Fˆ（X，1
T0

それが ∑bd3vcCab（Fa，Fb）＝－∇・Γac と表せることを使っている［10］．また上式の磁気面平均〈⋯〉を取ると，（2．15）

ˆˆ 0 /T0＋hˆ であり，それ以外が f0
この表式において f＝－eϕFˆに当たる．f の第一項，第二項はそれぞれ Fˆ 1 のポテンシャ∂∂（V ′〈nac〉）＋ （V ′〈（ΓaNC＋Γa＋Γac－nacus）
・∇s〉）＝0∂t∂s（2．10）

ル揺動 ϕˆ に対する断熱応答，非断熱応答部分を表している．f0 の中には ρ を通じてジャイロ角 ξ に依存する部分がを得る．ここで us ＝ ∂X/∂t は背景磁場の変化による磁気あるが，それは（2．13）式下の定義より O（δ1）までの（n，P）

面 s の移動を表している．また，磁気面内の体積を V としに寄与しない（速度空間の ξ- 積分で消える）ことがわかる．
307

Journal of Plasma and Fusion Research

特に圧力テンソルのアンサンブル平均については，
P ＝md vv fˉ ，P ＝md vv fˉ /2，

Vol.100, No.7

July 2024

場の形成過程における過渡的応答を扱う際に重要となる．
なお（2．17）式の反磁性フローに対し，∇・（u⊥1＋u‖1）＝ 0，

P＝m d3vvv（F0＋fˉ1）＝Ip0＋P‖1bb＋（I－bb）P⊥1，∇・（q⊥1 ＋ q‖1）＝ 0 を満たす 1 次のパラレルフローが存在する．このパラレルフローは以下のように表わすことが

‖1

3

2
‖ 1

⊥1

3

（2．16）

2⊥ 1

できる．
〈u‖1B0〉Xa1 ˜
u‖1，a＝   
B0＋  U，ea〈B02〉
〈q‖1B0〉5pa0 Xa2 ˜
B0＋    U，
q‖1，a＝   
2 ea
〈B02〉

T（s）
＝n（s）
は MHD 平衡 J0 × B0 ＝と書き表せる．p（s）000∇p0 を満たす磁気面関数のスカラー圧力である．一方，
1 次の非等方圧力はドリフト運動論方程式（2．4）式の解 fˉ1から決まるもので，Chew–Goldberger–Low 形式（CGLform）と呼ばれる磁力線方向と垂直方向の 2 成分から表さ

Xa1＝－

れる．（2．14）式第 1 項は ∇p0 からの寄与はなく，非等方′p0a－eaΦ′，Xa2＝－T0a′ ，
n0a（2．18a）（2．18b）

˜（2．18c）（ B ）＝B ×∇s・∇（ B ），〈BU〉＝0．

圧力テンソル ∇P1 によって駆動される．これが新古典輸

B0・∇送に対応する部分である．
（2．14）式第 2 項は E × B ドリフ

˜
U

0

1

0

2
0

˜ の具体的な求め方については［14］を参照せよ．また，U
˜に（2．18b）式 で ′＝d/ds を 意 味 す る．u‖1，q‖1 の う ち Uトによる流れであるが，E0＝－∇Φ0＝－∇〈ϕ〉ens は 2．2 節
で見たように Φ0 が磁気面関数であるとみなすと，この項比例する部分は MHD 平衡の理論にも出てくる Pfirsch-は Γ s に寄与しない．次に Coulomb 衝突項による摩擦項 R1

Schlüter フローに相当するものであるが，これは B0 を掛に目を向ける．衝突項を（2．5）式の線形化オペレータ CL でL
δ1 次の項まで考えるとすると，Maxwell 分布 F0 は CabにL寄与しない．また，摂動分布 Fˆ 1 について〈Cab（fˆ）〉ens＝0 でけて磁気面平均すると消える．磁気面平均部分はドリフ
ト運動論方程式を解くことによって求まる新古典輸送のfv）のうち－ F0 ρ・S0 ＋ fˉ1あるから，R1 に寄与するのは（x，の項のみである．ここで，fˉ1 はジャイロ平均された分布関寄与であるが，これは後で考える．これまでの議論をまとめて（2．14）式を整理すると以下のようになる．

L ˉf1）＝R‖1 と，磁力線平行方数であることから md3vvC（〈〉〈〉
K ×B ・∇s∂ n＋〈     〉＋〈 （ b×u）
・∇s〉
eB
∂t Ω

B0×∇s・∇・P1
R1⊥×B0・∇sΓs≡〈nu・∇s〉＝－       ＋      eB02eB02

向（パラレル方向）の摩擦力しか生まないことになる．従って fˉ1 は（2．14）式第 3 項に寄与しない．一方，有限 Larmor
半径効果の最低次の項である－ F0 ρ・S0 を CL に入れると

1



磁力線に垂直方向の摩擦力が生じ［13］，以下の関係式を

0
2
0

pl
＝ΓNC＋ΓC＋Γ ＋Γ⊥．

得る（簡単のため T0a ＝ T0b と置く）．

（2．19）

n0ama cR⊥1，ab＝－     bτab eB0

このように，新古典，古典，乱流の粒子束をジャイロ中∇p0a ∇p0b
3
1－zama/zbmb
×    －   －∇T0      zan0a zbn0b
1＋ma/mb
2za

ローとして［（2．14）-（2．19）式］表現することができた．

n0ama＝－   b× u⊥1，a－u⊥1，b
τab


－

（

（）

心のフローとして［（2．9），（2．10）式］，及び粒子位置のフ
3 種類のフラックスはここで見たように，ドリフト運動論方程式の解 fˉ1，背景プラズマの勾配，及びジャイロ運動論方程式の解 Fˆ 1 によって定義され，それぞれ明確に分離

）

ma
mb
3
q⊥1，a－
q⊥1，b     ，（2．17a）
mb
ma＋mb
5

されている．R1⊥ に比例する古典拡散や（2．17）式の反磁˜ に 比 例 す る Pfirsch-Schlüter フ
性 フ ロ ー，（2．18）式 の U

b×（∇p0a＋eana0∇Φ0）u⊥1，a＝           ，n0amaΩa5p0ab×∇T0a
q⊥1，a＝       ．2maΩaローは，MHD 平衡磁場と背景プラズマの勾配 Xa1，Xa2
のみで決定できる．これらは 2．1 節の輸送現象の分類に従えば流体的な輸送現象と言える．また，（2．19）式のフ（2．17b）

ラックスはいずれもドリフトオーダリングで δ2 次の項とみ出す粒子及び熱の流れで，diamagnetic（反磁性）フロー

なっているが，それを求めるのに必要なのは δ1 次までの分布関数（ fˉ1，Fˆ 1）であることは重要な点である．なお，こと呼ばれ，温度と圧力が磁気面関数であることから磁力

こでは粒子フラックスの表式にしぼって話を進めてきた

線に垂直かつ磁気面に沿って流れるため，径方向輸送には寄与しない．なお u⊥1，a の表式には摩擦力には寄与しな

mv2
が，（2．11）×    v のモーメントを取ることで径方向の2

い E×B フ ロ ー も 含 め て い る．（2．17）式 よ り，R⊥1，ab は

熱輸送方程式も同様に導かれ，熱フラックスについても反磁性フローに起因する垂直方向の摩擦力に対応してい（2．19）式と同様の分解ができる．

ここで，u⊥1，a，q⊥1，a は分布関数の－ F0 ρ・S0 の部分が生（ ）ることがわかる．（2．14）式の第 4 項は乱流場と分布関数

このように，新古典輸送と乱流輸送はそれぞれを記述

揺らぎの相関 K1 に駆動される，乱流輸送を表現している

する運動論方程式の上からも［（2．4），（2．7）式］，フラック部分である．最後に，（2．14）式の最終項は分極ドリフト

スの定義からも［（2．9），（2．19）式］明確に分類できること

による粒子束を表すもので，後に見るように両極性径電が示された．しかしながら，それは必ずしも 2 つの輸送308

Special Topic Article

2. Overview of Neoclassical Transport Theory	

過程が相互に影響しないということを意味しない．これS. Satake〈B0×∇s・∇・P1
ΓNC＝－       eB02までのほとんどの研究において両者を独立に取り扱ってきたのは，あくまでも分離できると仮定していたからで

〈

〉

（

）〉

B0・∇・P1
G 〈B0・∇・P1〉
B02＋      
＝－       1－  
22e〈B0 〉
eι
〈B02〉
〈B0 〉

あり，その理由は個々の輸送現象についての理解のためという面もあるが，実際上 2 つを分離しないと今までの計算機ではシミュレーションが不可能であったからとい  ＋

う側面も大きい． 2 つの輸送現象の相互作用のパスとして〈eζ・∇・P1〉eι＝ΓBP＋ΓPS＋Γna．

まず考えられるのは電場 E である．ヘリカル系の新古典輸（2．22）

送ではこの後で説明する両極性条件から大域的な径電場

ΓBP，ΓPS はそれぞれバナナ - プラトー，Pfirsch-Schlüter

E0＝－∇Φ（s）
が決定されるが，この径電場分布が乱流輸
0

フラックスと呼ばれるもので，それぞれ低衝突領域（衝送に影響を与えることが考えられる．一方，乱流の非線突周波数 ν ＜ R/vth）及び高衝突領域（ν ＞ R/vth）で主要形発展の結果，小半径方向に強いシア―を持つゾーナル
フローを生み出すポテンシャル構造〈ϕ˜ 1〉ens が発達する場の磁力線方向及びトロイダル方向成分の磁気面平均であ合があり，これがドリフト軌道を変化させて新古典輸送

り，そのような項が顕著になるのは圧力テンソル P1 に大な新古典輸送を生む．なぜなら，ΓBP と Γna は新古典粘性に影響を与える可能性もある．また，ジャイロ（ドリフト）

きな非等方性（P‖1－P⊥1）が生じる場合，即ち Coulomb 衝

運動論を解いた結果生じた磁気面上の密度や静電ポテン

突によって速度空間における捕捉 - 非捕捉軌道の境界近傍での分布関数の非等方性を Maxwell 分布に戻そうとする

シャルの偏りが，ドリフト（ジャイロ）運動論を解く際に高次の修正項として効いてくる場合も考えられる．こ効果が弱い場合である．一方 ΓPS は衝突頻度が高く荷電粒

うした相互作用の可能性を考え，（2．4）式と（2．7）式を分子が磁力線に沿ってトーラスの内側，外側を自由に往来

離せずにまとめて数値シミュレーションとして解く試み

できない場合に，磁力線方向に生じる圧力差によって駆が最近になってスーパーコンピュータの発達に伴い発展動される．また，Γna は新古典トロイダル粘性に比例する

してきている．その内容については 4 章で紹介する．

部分であり，その形から軸対称系ではゼロとなることがわかる．つまり，この項が非軸対称系の新古典輸送を特2．4

徴づける部分である．ここで ΓBP ＋ ΓPS の部分に（2．20）式

新古典輸送理論の詳細を使うと，次式が得られる．

以降，K1 の項を落として新古典輸送に話を限定し，もう少し理論の詳細について見ていくことにする．まず，

〉 〈（
1
1 ∂
＝－ 〈enE ＋R 〉＋ 〈 L 〉．
eι
eι ∂t

ΓBP＋ΓPS＝－（2．13）式と B0 の内積を取り，磁気面平均することで，δ1 次のパラレル方向のモーメンタムバランスの式が得られる．

〈

1‖ζ

〉

∂
〈B0・∇・P1〉＝e〈nE1‖B0〉＋〈R1‖B0〉－（mnu‖B0）．
∂t（2．20）〈）〉

G enE1‖＋R1‖ G ∂ mnu‖      ＋       
eιeι ∂t B0B01‖ζ（2．23）‖ζ

ただしここで G/B0＝b・eζ を使い，L‖ζ ＝ mnu‖・eζ 等とする．一般的な ΓBP，ΓPS の定義では E1‖ζ に関する項を分離し，また ∂L‖ζ/∂t の項は定常を仮定して無視することが多

なお，E1‖ は（2．13）式では E0 として現れていたものより高次のアンサンブル平均電場で，例えばトカマクの誘導電

いが後に考え直すことにし，∂L‖ζ/∂t の項を除いて各粒子場のような外部駆動電場を表すものとする．また，これ種の電荷 ea を掛けて和を取ると，

以降 B0 の時間変化は O（δ2）以上で無視できるものと仮定BPする．次に，（2．13）式と eζ ＝ ∂x/∂ζ の内積より，以下の関



〉〈

〈

1

a

係式が導かれる．

〈

PS〉＝－ ∑ eanaE1‖ζ＋R1‖ζ，a ＝0，
（2．

24）
a
a ＋Γa ）∑ e（Γιaとなり，磁気面を横切る電流を生まない（本質的両極性，〉∂∂Lζ ≡ （mnu・eζ）∂t∂t

intrinsic ambipolarity）ことがわかる．なお，ここで準中性条件 ∑aeana＝0 および Coulomb 衝突項のモーメンタ

pl＝eι（ΓNC＋Γ⊥）－〈eζ・∇・P1〉

ム保存性 ∑aR1‖ζ，a＝0 を使った．この性質は次に考える両＋e〈nE1‖ζ〉＋〈R1‖ζ〉．（2．21）

極性条件で重要となる．

この式はプラズマのトロイダル回転の時間発展を表し，

pl
さて，（2．19）式で分極フラックス Γ⊥を特に説明なく導トカマクにおいて誤差磁場等の影響による軸対称の破れ

入したが，これについて考察してみる．b ×（2．13）式よりに起因する新古典トロイダル粘性〈eζ・∇・P1〉がトロイダ

垂直方向フローの式を導くと，以下のようになる．

ル回転に与える影響の研究［12，15］などに使われる．なお（）

nE0×B0 B0×∇・P1 R1⊥×B0 ∂ nnu⊥＝   b×u⊥ ．
＋    
＋   
＋
B02
eB02
eB02
∂t Ω

（2．21）式では古典フラックスに関係する ΓC と R1⊥ の項は無視できるとして落としてある．次に，Boozer 座標系で

B0×∇s ＝（GB0－B2eζ）/ι と変形できることを使うと，ΓNC

ここで今，輸送オーダリングより最終項の時間微分の中を次のように分解できる．

で速く変動するのは E×B ドリフトの静電ポテンシャルの時間変化由来の寄与のみであると考える．つまり，これ309

Journal of Plasma and Fusion Research

を分極フラックスとして以下のように考える．
pl
nu⊥
≡

（

） （

July 2024

〈

〉

mana
∂Es
  ＝－∑ eaΓaNC＋   2 |∇s|2    ，
B
∂t
a
0

）

n
∂ n
∂Φ0' ∇s×b
b×u⊥ ≃ b× －      
B0
Ω
∂t Ω
∂t

mn ∂Φ0'＝－    ∇s．eB02 ∂tVol.100, No.7〈

∴ϵ0 |∇s|2＋（2．25）〉

c2∂Es|∇s|2   ＝－∑ eaΓaNC，
vA2
∂t
a

ϵ0c B022

vA2＝∑     ．a mana

pl
pl
＝〈nu⊥
・∇s〉
ここで Φ0' ＝ dΦ0/ds．これより，
（2．
19）
式の Γ⊥

（2．28a）（2．28b）に対応していることがわかる．一方，（2．23）式で現れたL‖ζ に対応して，L⊥ζ ＝ mnu⊥・eζ という量を導入すると，なお通常 c2/vA2≫1 であり，分極電流の効果により径電場の次のことがわかる．

される．（2．28a）＝ 0 を満たす径電場を両極性径電場と呼振幅は単純に eaΓaNC だけを考慮した場合よりかなり抑制

ぶ．非軸対称系プラズマでは両極性条件 Js ＝ 0 を満たす径mn ∂Φ0'∂ιmnΦ'L⊥ζ≃    B|∇s|2，
0×∇s・eζ＝－   
B02
B02 ∂t
∂t

〈

〉

∂
pl
L⊥ζ ＝eιΓ⊥．∂t

電場の解が， 1 つだけでなく複数個現れる場合がある．この問題については 3 章でその例が示される．一方，定常（2．26）

状態を仮定してドリフト運動論方程式を解く手法の場合∂
L ζ ≡ eιΓ‖pl と書くのがすると（2．23）式の最終項も，
∂t ‖
∂
∂
∂
pl
Lζ ＋
L⊥ζ ＝
Lζ ＝〈eιΓ‖pl ＋ Γ⊥
妥当であり，
〉
∂t ‖
∂t
∂t

当に振りながら ΓaNC（もしくは Γana）を求め，Js ＝ 0 を満

∴

〈

pl

〈

〉〈

〉

〉〈

pl
は Γ⊥
a ＝ 0 であるから，Es の値を入力パラメータとして適〉たす条件を探索するのが一般的な方法となる．トカマクなど本質的両極性が成り立つ配位では径電場は（2．28）式＝ eιΓ と書けることがわかる．つまり，径方向の分極フのように決めることはできず，径方向力学平衡（Radial

ラックスは磁気面平均したトロイダル方向の運動量の変force balance）と呼ばれるプラズマのトロイダル回転，化に比例して生じるものと解釈できる．これまでの結果

圧力勾配と径電場の関係式から決定される［17］．実際にを踏まえて，磁気面を横切る電流を求めてみると，以下は外部からのモーメンタムソース・シンクの情報が必要のように書けるとなる他，微視的乱流による粘性やモーメンタム輸送な
ど複合的な要素によって決まると考えられる（自発回転，

pl
pl
NC
na
Js≡∑ eaΓas＝∑ e（Γ＝∑ e（Γ
aa ＋Γ⊥a）a
a ＋Γa ）a＝aintrinsic rotation［18］）．

a

〈

〉

1
∂
∑ 〈eζ・∇・P1〉a＋ ∂t Lζ ．
ι a

さて，ここまで径方向のフラックスにのみ注目してきたが，新古典輸送理論のもう一つの重要な応用は温度密（2．27）

度勾配に駆動されるパラレルフロー，すなわちブートス
トラップ電流の評価である．その関係を含めた新古典輸

ここで，左から 3 番目の表現において（2．19）式にもともCとあった古典拡散と乱流輸送 Γ ＋ Γ の寄与が現れていな送の表式を導くために，今一度（2．20）式のパラレル方向

いが，これは新古典輸送だけ考えて無視しているわけでのモーメンタムバランスの式に戻ることにする．これ以はなく，この 2 つの項はそれぞれ，衝突項のモーメンタ∂
降， の項は無視できると仮定する．また，（2．20）式は∂t
mv2 5（2．13）
・mv‖b0 より得られたが，同様に mv‖   － b0
2T0
2



ムバランスおよび準中性条件［16］から径方向電流を生ま（ないため実際に寄与しない．同様に，（2．24）式でも見たように（2．27）式の左から 3 番目から 4 番目の表現に移るを内積して得られるモーメント量も合わせて考えると，際，ΓNC に含まれていた ΓBP ＋ ΓPS の項は本質的両極性の次のような関係式が得られる［19］．

ため消えている．このように，Js による両極性径電場の形成は新古典輸送現象による，しかも非軸対称系に特徴〈B・・P1〉－e〈nE1‖B0〉〈R1‖B0〉0 ∇            ＝    〈B0・∇・Θ1〉〈R2‖B0〉 aa的な現象であることがわかる．
それでは，Js から径電場 Es ＝－ Φ0' の時間発展の式を導〈Bu‖b〉
ab
ab
l11
－l12
＝∑            ，

2
ab
ab
b
l22
－l21
〈Bq‖b〉
5p0b
L ˉR1‖a＝ma ∑ d3vv‖ Cab（f1a），

く．この際，確かに径方向電流に実質的に寄与するのは

Γna と Γapl であるが，時間発展型のドリフト運動論方程式
pl∂（2．29）

5
R ＝m ∑ d vv x － C （fˉ ），（ 2）〈 ∂t L 〉/eι を高精度に評価するの'コードにおいて，Γ‖ ＝

）

b

‖ζ

∂Φ0
pl
は難しい．一方，Γ⊥a はそれ自体が  に比例する項とし∂t2‖aて定義されているので，磁気座標系における Poisson 方程

a

3

b

‖

2
a

m a v2xa2＝   ， P1a＝（2．16），
2T0a

式から，以下のようにして径電場の時間発展を求めるの

 （が一般的である．

（）

v2
5
Θ1a＝ma d3v v‖2－ ⊥ xa2－
2
2

∂Es
pl
NC
ϵ0〈|∇s|2〉   ＝－∑ e（Γaa ＋Γ⊥a）∂taLab（2．30）

1a 31）（）b b － 3 I）fˉ ．（2．
0 0

1

1a

なお（2．
29）式はジャイロ中心分布関数 fˉ1a を Legendre 及びLaguerre 多項式で展開したうちの低次の項のみで近似し310

Special Topic Article

2. Overview of Neoclassical Transport Theory	S. Satakeて線形化衝突項に入れた場合に得られる関係（13 モーメ
ント近似［20］）である．（2．29）式は磁気面平均した磁力式となっている．さて，（2．
22）
式から〈B0・∇・P1〉∝ ΓBP〈Bu‖a〉/〈B2〉Na1 Na2Γbn＝ΓBP＋Γna      
＝                2qbn＝qBP＋qna a
Na2 Na3   〈Bq‖a〉〈B2〉5p0aが 示 さ れ て お り， 径 方 向 の 新 古 典 熱 フ ラ ッ ク ス を〈d v （x － 2 ）fˉ v ・∇s〉と 定 義 す る と，

qNC ≡ T

3

5

2

c

Na1 Na2 Xa1＋       
．（2．33a）Na2 Na3 Xa2



線方向の粘性と摩擦力，フローの関係を結び付ける関係

1 d

（2．
21）
式－
（2．
23）
式と同様の計算からqNC ＝qBP ＋qna ＋qPS，

La1 La2 Xa1＋       
．（2．33b）
La2 La3 Xa2



qBP∝〈B0・∇・Θ1〉も示すことができる．さらにここで径方向フラックスについて，バナナ－プラトーと非軸対称成分の和と，Pfirsch-Schlüter の 2 つに分離することを

Nai，
Lai）
は ˆ a の解（より正確には ˆ a を
ここで係数行列（Mai，考える．すなわち，ΓNC＝（ΓBP ＋ Γna）＋ΓPS ＝ Γbn ＋ ΓPS，

Xa1，Xa2，〈Bu‖a〉，〈Bq‖a〉の線形和で表した時に
（2．
32a）
式

NC

BP

na

PS

bn

PS

PS

＋ q ）＋q ＝ q ＋ q ．説明は省くが，Γ ，
˜ 1‖〉
˜ 2‖〉
〈UR
〈UR
qPS は ΓPS ＝－    ，qPS ＝－    と表せるため［19］，eeを満たすように決まる係数）に関する速度空間積分及び磁これも（2．29）式同様にパラレルフローと関係づけることに対する式として整理すると，以下のような方程式が得ができる．そうすると，（2．29）式の関係式を使って粘性，られる．
（簡単のためここではイオン種は 1 種類としたが，フロー，径方向輸送と駆動力（（2．18b）式の Xa1，Xa2 およ

複数イオン種の場合も同様の表式が得られる）q＝（q

気面平均から決定される．
（2．
29）
式と
（2．
33）
式を結び付け，
Γbn，qbn，JBS ≡〈J‖B〉/〈B2〉＝〈（einiu‖i － eneu‖e）B〉/〈B2〉び XE ≡〈E‖B〉/〈B 〉）を結びつける式を導けそうである．
2

Γebn

ee
ee
ei
ei
e
L11 L12 L11 L12 L1E
bn
ee
ee
ei
ei
e
qe
L21 L22 L21 L22 L2E
ie
ie
ii
ii
i
bn
  
Γi ＝            L11 L12 L11 L12 L1E
ie
ie
ii
ii
i
L21 L22 L21 L22 L2E
qibn
e
e
i
i
LE1 LE2 LE1 LE2 LEE
JBS

ここでは，モーメント法と呼ばれる手法の概略だけを示す．まず，出発点はドリフト運動論方程式（2．4）式を独立
変数を（x，v，ξ ＝ v/v‖）として書き直した以下の式である
［19］．

（

）

5
⎧
⎫
（V‖－CaPAS）ˆ a＝vda・∇s⎨Xa1＋ xa2－ Xa2⎬
⎩
⎭
2


＋

fˉ1a
ea
a≡  －
Fa0
Ta

  l



dl
（E‖B－B2XE），
B

ˆ a＝－（l＝1）．

Xe2
Xi1 ．

（2．34）Xi2
XE

ab

係数（Mai，Nai，Lij ）を具体的に求める方法としては主に変分法を使う方法［21-23］と Monte-Carlo 法を使う方法

ea
XEBvξ－（V‖－CaL）（l＝1），（2．32a）
Tav（1－ξ2）∂V‖＝vξb・∇－    b・∇lnξ
，
2
∂ξ

Xe1

［24］が知られているが，ここではその詳細には触れない．
（2．34）式のように表現すると，新古典輸送（径方向粒子，
熱フラックス，ブートストラップ電流）が（Xa1，Xa2，XE）（2．32b）

を駆動力として引き起こされる輸送過程であることがよくわかる．特に，粒子フラックスが Γbn が圧力勾配と径電（2．32c）

場 Xa1 だけでなく温度勾配 Xa2 や誘導電場 XE に駆動されたり，熱フラックス qbn が Xa1 にも駆動されるといった輸送（2．32d）

ab

行列 Lij に非対角成分があること，また電子（イオン）が自身の勾配だけでなく摩擦力を介してイオン（電子）のこ こ で，（2．32c）
式の積分は磁力線に沿った積分であり，（l）は  を Legendre 多 項 式 P（ξ）
［P0 ＝ 1，P1 ＝ ξ，l

勾配にも駆動されることなど，新古典輸送の重要な性質

P2 ＝（3ξ － 1）/2，…］で展開ししたときの l 次の成分を意味

がこの 1 つの式にまとめられている．一方で（2．22）式，する．また，左辺では衝突項をピッチ角散乱オペレーター（2．23）式のように表現した場合，新古典輸送は圧力の非

2

a

νD ∂
∂
a
CaPAS ＝   （1 － ξ2） で近似している（νD は 90 °散乱
2 ∂ξ∂ξ

等方性（＝粘性），あるいはパラレル方向の摩擦力の結果

衝突周波数）．衝突項に関するこの近似は一見大胆に見え（2．20）式と（2．29）式に現れるように Coulomb 衝突による生じると見ることもできる．いずれの見方においても，
パラレル方向の摩擦力，モーメンタムバランスが新古典L ˆ a）＝ d3v mv‖CPAS（ˆ a）＝0というモーるが， d3v mv‖C（

輸送の評価に重要である．なお，Pfirsch-Schlüter フラックスについても同様に表せ［25，26］，以下の形となる．

メンタム保存性を数学的に満たすことが知られている．
また，（2．29）式において左辺のテンソル量 P1，Θ1 は （l ＝ 2）ababl11 －l12 Xb1
ΓaPS
1
˜2
  
．
ab ab
PS ＝－〈U 〉∑           
e
e
qa
Xb2
b a b －l21 l22

成分で決まる一方，右辺のパラレルフロー u‖，q‖ は （l ＝ 1）

（2．35）

成分のみで決まるという関係になっている．このような関係を使って a の解を求めると，パラレル粘性，フロー，ab˜2〉は MHD 平衡磁場から決まり，摩擦係数 ljk
も線形化〈U

径方向輸送，勾配 Xa1，Xa2 を結び付ける線形の関係式が衝突項の性質のみから決定されるため，（2．35）式はドリ得られる［19］．

フト運動論方程式を解かずとも決定できる．そのような
意 味 に お い て，Pfirsch-Schlüter フ ラ ッ ク ス は 2．1 節 の〈Bu‖a〉/〈B2〉Ma1 Ma2
〈B・
・P1〉
0 ∇
      ＝                
2
〈B・
∇
・
Θ
〉
Ma2 Ma3   〈Bq‖a〉〈B2〉
0
1
a
5p0a

分類上流体的な輸送過程と見ることができる．Pfirsch311

Journal of Plasma and Fusion Research

Vol.100, No.7

July 2024

Schlüter フラックスは磁気面上のパラレルフローの変動
˜
部分（∝U）によるパラレル方向の摩擦力の変動（＝パラモーメント法の説明の中で衝突項の性質については特に

B02
レル粘性の変動，（2．22）式の 1 －  に比例する項）に起〈B02〉

ljk は自己随伴性から ljk ＝lkj ，またモーメンタム保存性から ∑al1k ＝0 を満たす．後者の性質が ΓBP と ΓPS の本質的両因する径方向輸送であるため，Pfirsch-Schlüter フローと極性を保証するが，自己随伴性はここでは特に必要とは同様に流体的な輸送過程になっているのはつじつまが合されない．しかし，新古典輸送における Onsager 対称性

うし，Pfirsch-Schlüter フラックスは衝突周波数が高い場というもう一つの重要な性質は，この自己随伴性から導合に顕著になるから，その観点からも流体的な描像と合
˜2〉は磁場のトーラス性を表す因子になっう．もちろん〈U

かれる．Onsager 対称性は（2．34）式において，触れなかったが，（2．29）式，（2．35）式において摩擦係数
ab

ab

ba

ab

ab

ており，故に ΓaPS と qaPS は新古典輸送の一部である．

ba

ab

ba

a

a

Ljk ＝Lkj ，Lj3 ＝－L3j ，LjE＝LEj（ j，k＝{1，2}）（
 2．37）

最後に，トカマクや新古典最適化配位の一種である準対称配位において重要な性質である本質的両極性（intrinsicの形で現れる．この性質は新古典輸送がエントロピーを増

ambipolarity）の概念について説明する．既に（2．24）式

大させる方向にのみ起こる散逸過程である事（Boltzmannで見たように，衝突項のモーメンタム保存性および準中の H- 定理）とも関係しており，文献［25］によくまとめら性条件より，径方向の新古典粒子フラックスのうち，パれている．重要な補足として，古典輸送と乱流輸送も同

ラレル粘性に関係する ΓBP と ΓPS については径方向の電様に Onsager 対称性を満たす輸送行列によってフラック

流を生まない．よってトロイダル粘性由来のフラックス

スと駆動力が結び付けられる［30］．なお，トカマクにおΓnq∝〈eζ・∇・P1〉及び分極フラックス Γpl だけが径電場の時間発展に関与することが（2．27）式で示された．ここ

いては新古典フラックスが本質的両極性を満たすだけでなく，径電場－ Φ0' に依存しなくなることも重要な性質ででは eζ として Boozer 座標系における任意のトロイダル角

あるが，最近の論文［31］において，この性質を計算で再を取っているが，新古典粘性をパラレル方向ともう 1 つ

現するには衝突項のモーメンタム保存性だけでなく自己の方向に分離する方法は 1 つではなく，ポロイダル成分，

随伴性も必要であることがわかりやすく示されている．

あるいは任意の成分を取ることも可能である［12，27］．もし，磁場が準軸対称，準ポロイダル対称，あるいは準ヘ

以上，本章ではドリフト運動論方程式とジャイロ運動

リカル対称などの準対称配位であれば，磁場に平行な方論方程式の導出から始めて，新古典輸送理論の基礎につ向と対称性の方向（＝ Q とする）の 2 成分に新古典粘性

いて概観してきた．次章では，具体的に新古典輸送の計を分離した場合，〈Q・∇・P1〉＝ 0 となるのでトカマクと算を行う上でどのような仮定や近似が使われているのか，BP

同様に新古典輸送は本質的両極性を満たす．一方，Γ

それが時代とともにどう進化してきたかについて見ていと

ΓPS（さらに ΓC も）の磁場配位に依らない本質的両極性くことにする．

は，Coulomb 衝突項のモーメンタム保存性から導かれた．
ここで線形化 Coulomb 衝突項，特に Landau の Fokker-

参考文献

Planck オペレータ［28，29］の持つ重要な性質を述べておこLう．（2．5）式の Cab（ fˉ1a）は，［ 1 ］ITER Physics Basis Chap. 2, Nucl. Fusion 39, 2175
(1999).
［ 2 ］Y. Yamada et al., Nucl. Fusion 45, 1684 (2005).
［ 3 ］A. Ishizawa et al., Nucl. Fusion 57, 066010 (2017).
［ 4 ］F. Warmer et al., Phys. Rev. Lett. 127, 225001(2021).
［ 5 ］R. Sakai et al., Fus. Eng. and Des. 149, 111322(2019).
［ 6 ］A.J. Brizard and T. S. Hahm, Rev. Mod. Phys. 79, 421(2007).
［ 7 ］H. Sugama, Phys. Plasmas 7, 466 (2000).
［ 8 ］プラズマ・核融合学会 編：プラズマシミュレーション（京都大学学術出版会，2018）第 2 章．
［ 9 ］A.H. Boozer, Phys. Fluids 24, 1999 (1981).
［10］H. Sugama, Rev. Mod. Plasma Phys. 1, 9 (2017).
［11］H. Sugama and W. Horton, Phys. Plasmas 2, 2989
(1995).
［12］S. Satake et al., Plasma Phys. Control. Fusion. 53,
054018 (2011).
［13］F.L. Hinton, Basic Plasma Physics 2 (edited by
A.A. Galeev et al., Handbook of Plasma Physics)(North-Holland Publishing Company 1984) Chap.
1.5.
［14］D.A. Spong, Phys. Plasmas 12, 056114 (2005).

L ˉ
Cab（fa1）＝Cab（ fˉa1，Fb0）＋Cab（Fa0，fˉb1）
TF＝Cab（ fˉa1）＋Cab（ fˉb1）
Tとフィールド粒子衝突項のようにテスト粒子衝突項 Cab

FCabに分けられる．これらは以下のような関係式を満たす．

・粒子数保存則：d vC （ fˉ ）＝d vC （ fˉ ）＝0，
T
ab

3

a1

3

F
ab

b1

（2．36a）

・モーメンタム及びエネルギー保存則：





T
F
ma d3v{v，v2}Cab（ fˉa1）＋mb d3v{v，v2}Cba（ fˉa1）＝0，（2．36b）

・自己随伴性（Self-adjointness）：
ˉ 1 T ˉfˉ1 T
d3v Cab（ f1）＝ d3v  Cab（ˉ 1），
Fa0
Fa0



Ta



fˉa1 F ˉfˉb1 F ˉ
d3v  Cab（ fb1）＝Tb d3v  Cba
（ fa1）．（2．36c）
Fa0
Fa0





なおここで {v，v2} は v または v2 のいずれかを意味する。
312

Special Topic Article

2. Overview of Neoclassical Transport Theory	S. Satake(2009).
［25］H. Sugama and W. Horton, Phys. Plasmas 3, 304(1996).
［26］J.D. Lore, Measurement and Transport Modelingwith Momentum Conservation of an Electron Internal
Transport Barrier in HSX (Doctor thesis, Univ. Wisconsin-Madison, 2010) Chap. 3.
［27］K.C. Shaing et al., Nucl. Fusion 50, 125012 (2010).
［28］L. Landau, Phys. Z. Sowjetunion (USSR) 10, 154(1936).
［29］M.N. Rosenbluth et al., Phys. Rev. 107, 1 (1957).
［30］H. Sugama and W. Horton, Phys. Plasmas 4, 405(1997).
［31］M. Honda, Phys. Plasmas 30, 092305 (2023).

［15］M. Honda et al., Nucl. Fusion 54, 114005 (2014).
［16］Sugama et al., Phys. Plasmas 3, 2379 (1996).
［17］J.D. Callen et al., Phys. Plasmas 17, 056113 (2010).
［18］J.E. Rice, Plasma Phys. Control. Fusion 58, 083001(2016).
［19］H. Sugama and S. Nishimura, Phys. Plasmas 9, 4637
(2002).
［20］R. Balescu, Transport Processes in Plasmas Vol. 2(North-Holland, Amsterdam, The Netherland, 1988).
［21］S.P. Hirshman et al., Phys. Fluids 29, 2951 (1986).
［22］W.I. van Rij and S. P. Hirshman, Phys. Fluids B 1,
563 (1989).
［23］H. Sugama and S. Nishimura, Phys. Plasmas 15,042502 (2008).
［24］A. Matsuyama et al., Phys. Plasmas 16, 052501

313

J. Plasma Fusion Res. Vol.100, No.7 (2024) 314-320

小特集 古くて新しい新古典輸送理論の新展開

3 ．新古典輸送計算の実際 ～近似モデルと新展開～
3. Neoclassical Transport Calculations in Practice:
Approximate Models and New Developments

3．1

ドリフト運動論方程式の数値的解法と種々の近似

3.1 Numerical Schemes for Drift-Kinetic Equation and
Several Approximations in them
佐竹真介

SATAKE Shinsuke
核融合科学研究所

構造形成・持続性ユニット（原稿受付：2024年 2 月21日）

2 章で導出したドリフト運動論方程式を実際に数値的に解いて新古典輸送を評価する際には，計算量を減らすために様々な近似や仮定が使われることが多い．ここでは，代表的な新古典輸送シミュレーションコードにおいて使われている衝突項やドリフト軌道などに関する近似について説明し，それらが計算結果にどのような影響を与えるかを説明する．
Keywords:
neoclassical transport theory, drift-kinetic equation, gyrokinetic theory, numerical simulation

3．1

ド リフト運動論方程式の数値的解法と種々
の近似

（

）

（

ma 3/2
1 maU 2exp －Fa0＝n0a      ＋μB0
T0a
2πT0a
2

2 章で見てきたように，新古典輸送ではドリフト運動論
方程式を解くことでガイディングセンター分布関数の局



所 Maxwell 分布からのずれを求め，それを用いて新古典）

eaΦ1・exp －   ．（3．1．1c）
T0a

ここで新たに入ったのは時間微分の項 ∂fˉa1/∂t と磁気面

粘性や新古典熱・粒子フラックス，ブートストラップ電流などを評価する．しかし，その計算を具体的にシミュ

上 の ポ テ ン シ ャ ル 変 位 Φ（s，θ，ζ）， 衝 突 項 の 非 線 形 項
1

レーションコードに実装する上では様々な近似や仮定がCab（fˉa1，fˉb1），そして vda・∇fˉa1 の項である．これらはドリ用いられる．本章では過去から現在まで使われている様々

フトオーダリングでは δ2 次の項という扱いになるが，あな新古典輸送計算法において，どのような近似が用いらる条件下では無視しえない影響を与えることもあるためれているかを解説し，それがどのような影響として計算残してある．また，（2．4）式では ∇Φ1 の項は E‖ の中に含結果に現れるのかについて見ていく．

まれていたが，ここでは Fa0 を（3．1．1c）式のように置くこ

ドリフト運動論方程式（2．4）式において，2．3 節では落とで E‖ からは分離されていることに注意．以下，それぞとしていた高次の項も含めて改めて書いてみる［1，2］．

れの項の扱いについて順にみていこう．
1 ）∂fˉa1/∂t の扱い：この項を含めて解くか，解かないか

∂fˉa1
1
∂
   ＋ Ub・∇－ b・（μ∇B0＋e∇Φ1）  ＋vda・∇ fˉa1
ma
∂t
∂U

（

はドリフト運動論方程式（3．1．1a）式を定常解に向かう時
間発展の問題として解くか，定常解を満たす fˉa1 を探索す

）

eaE‖U
L ˉ
＝－ vda・∇－     Fa0＋∑ Cab（fa1）＋Cab（fˉa1，fˉb1）,
T0ab
（3．1．1a）る問題として解くか，という数値解法の設定の違いに対応する．時間発展を解く場合，（2．28）式を用いて径電場の時間発展も同時に解くことが多い．また従って本質的両1
vda＝   b×［maU 2b・∇b＋μ∇B0＋ea∇（Φ0＋Φ1）］，
eaB0
（3．1．1b）極性（intrinsic ambipolarity）の成り立たない非軸対称系においてよくその立場が取られる．計算手法としては，
時間発展を解くコードでは分布関数を多数個の仮想粒子
̇ U）
̇ に（マーカー）をガイディングセンター運動方程式（X，National Institute for Fusion Science, Toki, GIFU 509-5292, Japanauthor’s e-mail: satake.shinsuke@nifs.ac.jp

314Ⓒ2024 The Japan Society of Plasma
Science and Nuclear Fusion Research

Special Topic Article

3.1 Numerical Schemes for Drift-Kinetic Equation and Several Approximations in them	S. Satake

従って動かし，そのマーカーの分布と重みから fˉa1 を表現しており，径方向のドリフトは無限小の仮想変位としてする粒子コードとして実装されることが多い．一方，定扱われることになる．このようなモデルを局所近似（local

常解を探索するコードでは径電場は入力値として与えられ，それを満たす解 fˉa1 を何らかの多項式展開，例えば 2．4

approximation）とも呼ぶ．コンピュータ資源が乏しかった時代に作られた新古典輸送計算法が基本的に局所近似

節で出てきた Legendre と Laguerre 多項式（速度空間）モデルであったのは，単一磁気面の計算に問題を落とし込

及び Fourier 多項式（配位空間）などで展開し，それらのめるという大きな利点があったからである．しかしなが係数が満たす関係を（3．1．1）式から導き，反復法や変分法

ら，実際の閉じ込め磁場中の荷電粒子のドリフト幅は背などで答えを求める，というものである．粒子法は 3 次元

景温度密度勾配と比べて無視できるほど狭いとは言えな配位への適用も容易で，数値アルゴリズムもドリフト運

い状況がよくあり（例えばトカマクのペデスタルやヘリ動論方程式の形から直接的に決まる一方，有限のマーカー

カルリップルに捕捉されたイオンの軌道など），この効果

で分布関数を表現することの不連続性や後に述べる衝突を取り入れた計算が重要となる．有限軌道幅効果を含め項をモンテカルロ法で実装することにより，数値ノイズ

た大域的新古典シミュレーションはスーパーコンピューが乗りやすくその抑制に数値計算上のテクニックが要求タの発達とともにこの 15 年ほどでようやく発展してきた．

,制される（例えば重み平均（weight-averaging）法［3］径方向ドリフトを含めた複雑な軌道を取り扱うために基御変量（control-variate）法［4］など）．多項式展開の場合，本的に粒子ベースのモンテカルロコードであることが多計算コードは一般に巨大な係数行列を扱う形になり，そいが，最近ではジャイロ運動論シミュレーションコード

の解法の収束性や安定性を得るための技法が必要となる．

が大域的シミュレーションとして発展する中で新古典輸2 ）磁気面上のポテンシャル変位 Φ（s，θ，ζ）：この効果1

送計算も含めて計算できるようになってきたため，ジャは従来の径電場しか考えない新古典輸送理論では無視さ

イロ運動論コードで使われる有限要素法など，分布関数

れていたが，近年高 Z イオンの不純物輸送を考える上で重を 5 次元位相空間内のメッシュ点での値を解く手法に基要視されてきた効果である．なお，ここで考える Φ1 は微づいた大域的新古典輸送計算コードも増えてきている．

視的乱流輸送を引き起こす〈ψˆa〉ξ（（2．7）式）とは異なり，また，vd・∇ の中には径方向成分（vd・∇s）だけでなく磁時間的にゆっくり（捕捉軌道のイオンのバウンス周期よ

気面上のトロイダル，ポロイダル方向成分も含まれる．ヘり十分遅く）変化する準定常的な成分で，イオンと電子

リカル系の局所近似新古典輸送コードの中では，この磁気のドリフト運動論方程式を解くことで現れる磁気面上のd ≡ vd－
面接線成分（v（vd・∇s）es）だけを残す近似を用い僅かな電荷の偏り ∑a ＝ i,e ean1a によって生じる磁気面上のるものが主流である．こうしたモデルは径方向ドリフト

静電ポテンシャルである．ただしこの場合も準中性条件の有限性だけを無視しているので，径方向局所（radiallylocal）モデルと呼ばれることもある．v
d の中には磁気ド∑a ＝ i,e〈ean1a〉＝ 0 が課される．Φ1 は eΦ1/T ～ O（δ）の微小量と想定され，実際の計算でも高々 eΦ1/T ～0.1 程度である．

リフト（grad-B と曲率ドリフト）による捕捉粒子の歳差

これが不純物イオンに対して重要視されるのは，ガイディ（precession）ドリフトや，E×B ドリフトの効果が入る．

ングセンターの運動方程式において

特に本質的両極性が成り立たないヘリカル系に特徴的な，
径電場に対する新古典輸送の強い依存性を調べるために，

vE1×B
ZeΦ1 R Ze∇‖Φ1ZeΦ1   ～   
，    ～   
vB×∇B
T a
TΔB
μ∇‖B

E×B ドリフト項のみ残した径方向局所モデルが存在し，
代 表 的 な も の と し て は DKES-PENTA［7-9］，GSRAKEと見積もられ（ΔB は磁力線に沿った磁場の強弱の振幅，［10］などが知られる．磁気ドリフトの磁気面接線成分もR 及び a はプラズマの大半径，小半径），電荷 Z が大きい含めた径方向局所モデルは比較的最近考案され，局所近ほど無視できなくなるためである［1］．Φ1 ポテンシャル

似版 FORTEC-3D［11，12］，SFINCS，KNOSOS［13］などを含めた計算をする場合，ジャイロ運動論でもよく使わが該当する．また，軸対称系に限定されるが局所近似でれるように断熱応答を仮定し， 0 次分布関数 Fa0 分布を（2．33）式の関係式に対する数値解法の計算コード（逆行（3．1c）式と置いた場合の fˉa1 に対するドリフト運動論方程列法 , Matrix Inversion［14］）などもこの類の計算法の一式を解く．その代表的な計算コードとしては EUTERPE

つである．磁気ドリフト項を含めた計算法（ただし下で［1］，FORTEC-3D［2］，SFINCS［5］，などが知られてお述べるバウンス平均近似を用いない場合）ではドリフト

り，いずれも 3 次元配位に適用可能である．また，NBI

運動論方程式における位相空間の保体積性（Liouville の加熱による高速イオンが作る Φ1 ポテンシャルに着目した定理）やエネルギー保存則が成り立たないため人為的なシミュレーションもある（GNET［6］）．3．3 節では Φ1 ポテ

ソース・シンク項が必要になるという問題点があるが［11］，ンシャルの不純物の新古典輸送への影響を調べた研究例

径電場ゼロ近くでの新古典輸送の径電場依存性がより正を紹介する．

しく表される．なお，上記の問題を解消できる新しい径3 ）非局所効果 vda・∇fˉa1：この項を含めた計算は，ガ方向局所モデルも提案されているが［15］，複雑な計算がイディングセンター軌道の磁力線を横切るドリフトの有必要となるためまだコードとして実装されたものはない．

限性を考慮して解くことになるため，非局所効果とか，また，ある限定的な状況にのみ妥当な近似を採用はして有限軌道幅効果などと呼ばれることがある．一方，このいるが，局所近似モデルに有限軌道幅効果を含めた計算項を無視したドリフト運動論方程式は磁気面ごとに独立モデルを構築する試みも近年なされている［16，17］．
315

Journal of Plasma and Fusion Research

Vol.100, No.7

July 2024

一 方 E×B ド リ フ ト 項 の み を 扱 う 場 合，vE×B ＝ E×

磁場配位に依らずブートストラップ電流を計算したい場

B/B → E×B/〈B 〉という近似をすることで位相空間の保合には，衝突項のモーメンタム保存性が必須である．ピッ

体積性が保たれるため［11］，この近似は非軸対称系の局

チ角散乱オペレータでは磁力線に平行方向のモーメンタ

所近似新古典計算法として従来からよく採用されている．

ムが保存しないため，補正を考える必要がある．一つの2

2T方法としては，テスト粒子衝突項 Cabとフィールド粒子衝

しかし，この近似は本来圧縮性流れである E×B フローにFT
突項 Cabという分離において，Cabをピッチ角散乱で近似

対する非圧縮近似に相当するため，E×B フローが大きい場合（ポロイダルマッハ数 Mp＝E/（Bpvth）が 1 以上）に起

F
T
し，Cab
として Cabで生じた変化をキャンセルする項を導こる新古典粘性，新古典輸送の径電場に対する非線形依

入するというやり方がある．これは初期の単一エネルギー

存性［18，19］を正しく扱えないことに注意が必要である．

近似の計算法でよく採られた方法である．また， 3 ）でまた，局所近似のドリフト運動論方程式の解法に関して名前を挙げた DKES コードは，ピッチ角散乱オペレータ

は，（θ，ζ）平面でのドリフト軌道をそのまま扱うものと，を使って（2．32a）式のような形式の，非圧縮 E×B フロー

磁場の弱い領域に捕捉された粒子軌道についてそのバウを入れた単一エネルギーのドリフト運動論方程式を解く．
この際，（2．32a）式右辺の（V‖ － CaL）（l ＝ 1） は無視されてンス平均したドリフト運動を扱うモデルがある．前者はDKES，後者は GSRAKE や KNOSOS，PENT［20］などがいるため，DKES が求める新古典拡散係数（エネルギー

代表的である．バウンス平均コードは他の計算法に比べモーメントを取ると（2．34）式の Lij になる）はこのままでて高速に新古典輸送計算を解くことができるが，磁力線は軸対称系において本質的両極性を満たさない．しかし，に平行方向の物理が入らないためブートストラップ電流（2．33）式に現れる係数（Mai，Nai，Lai）は，DKES の計算の計算はできない．また，捕捉軌道が衝突散乱される前結果を基にモーメンタムバランスを満たす形に修正してに十分にバウンス運動できることが定式化の前提となっ求めることができ，その関係を利用してピッチ角散乱近ているので，Pfirsch-Schlüter 領域の計算は別途行う必要

似を使った計算から軸対称系にも適用できる新古典輸送がある．なお，数値シミュレーションによる新古典計算

拡散係数を評価できる［22，23］．この関係を利用した新古が発展する以前からヘリカル系の新古典輸送を衝突周波

典輸送コードが，PENTA コード［9］である．

ab数，径電場の大きさ毎に区切って解析的に扱う計算式が衝突項の保存則を満たすより厳密なオペレータの実

多く考案されたが［21］，こうした解析モデルは基本的に

装 法 と し て， テ ス ト 粒 子 項 と し て，（2．36）式 に 挙 げ た

局所近似であり，捕捉粒子のバウンス平均した径方向ド

Landau オペレータのピッチ角散乱とエネルギー散乱を

リフト速度と，衝突や E×B ドリフトによる捕捉 ‐ 非捕捉

T
数値的に模擬するものを用意する．Cab
は速度空間におけ軌道の遷移の時間スケールから粒子拡散係数を評価する

る 2 階の偏微分方程式の形をしているので，数値的には多モデルになっている．

項式展開で表現された分布関数の微分系，差分，あるいはT
粒子シミュレーションでは Cabの Fokker-Planck 的な性質

4 ） 衝 突 項 に 対 す る 近 似：（3．1．1a）式 右 辺 に 現 れ る
Coulomb 衝突項の扱いも，新古典輸送シミュレーション法にを利用した乱数による速度空間でのランダムウォークとよって大きく扱いが分かれる部分である．まず，最も単純化して実装する［24］こともできる．これに対し，フィールされたモデルとして，衝突項を速度のピッチ角変数 ξ ＝ v‖/v

ド粒子衝突項は微積分方程式［25］の形をとり，数値的なνDa ∂∂（1 － ξ2） でに対する散乱オペレーター，CaPAS ＝
2 ∂ξ∂ξ

扱いが複雑になる．そこで，（2．36）式の保存則や自己随近似する方法がある．これは（2．32）式でも述べた方法でF
アドホックに Cabを決定するという方法が新古典シミュある．一般に本質的両極性が成り立たないヘリカル系でレーション，乱流シミュレーションの両方でよく採用さは，衝突項のモーメンタム保存性は小半径方向の輸送をれる［26，27］．特に最近では，多粒子種プラズマに対応し計算する上では重要ではなく，磁場リップルに捕捉されて Ta≠Tb の時に（2．36c）式の自己随伴性が満たされないとT
伴性を満たすように，数値計算された Cabの結果に対してた軌道と非捕捉軌道のピッチ角散乱による遷移が表現でFいう従来のアドホックな Cabの問題を解決したモデル［28，

きれば十分な精度で径方向輸送を計算することができる．

29］が多くのドリフト運動論コードとジャイロ運動論コー

また，衝突によるエネルギーの変化を無視することで，ドに実装されている．2．4 節で説明したとおり，線形化衝計算をより単純化できる．すなわち，径方向局所モデル

突項の自己随伴性は Boltzmann の H- 定理や，トカマクにでポテンシャルとして Φ（s）
のみ考えれば，ある磁気面に0

おいて新古典フラックスが径電場に依存しなくなる性質

束縛されて運動するガイディングセンターに沿っての運を数値的に再現する上で重要な要素である．しかしなが動エネルギーは Coulomb 衝突以外では変化しないため，ら，このようなアドホックなフィールド粒子衝突項には

ピッチ角散乱近似を用いることで，速度 v を変数ではなくパラメータとして扱うことができる．すると分布関数 fˉ1aやはり近似が入ってしまう．特に問題になるのは衝突周波数が高く Pfirsch-Schlüter フラックスが重要になる場の変数としては（ζ，θ，ξ）の 3 変数に落とすことができるたF
合で，アドホックな Cab
では新古典シミュレーションの定め，ヘリカル系でも必要な計算量を大きく減らすことが量的正確性が落ちることが指摘されている［30］．そこで，できる．こうした近似は単一エネルギー（mono-energy）

F
計算のコストはかかるが数値的に Cab
を元の定義式に従っ近似と呼ばれる．

て計算する，完全線形（full-linearized）オペレータを用一方，トカマクのような本質的両極性が重要な場合や，いた計算法も発達している［30，31］．フィールド粒子衝突316

Special Topic Article

3.1 Numerical Schemes for Drift-Kinetic Equation and Several Approximations in them	S. Satake

項の正確な扱いは，今後核融合炉内のタングステンの輸るには Js ＝ 0 が必要であるから，準対称配位では O（δ0）

送問題など，Pfirsch-Schlüter 領域の衝突周波数を持つとの速いフローを考えると矛盾が生じる．したがって，速想定される高 Z イオンを扱う計算において重要となる点でいフローが起き得るのは軸対称トカマクの場合に限られFある．アドホックな Cab
オペレータに対しても高衝突領域

ると予想される．さて，このような速い流れがある場合，での修正項が考案されており［32］，ジャイロ運動論コー

ドリフト運動論方程式は u0 に乗った座標系で記述されるドへの実装も進められている［33］．さらに，トカマクの［36］．この時の顕著な効果の一つとして，重イオン種のペデスタル部など急峻な背景の温度勾配がある状況では，密度がトロイダル回転による遠心力によって磁気面の外線形化衝突項の近似では不十分と考えられる．そのよう側に偏ることが知られている．即ち，ゼロ次の密度分布な状況に対応すべく，（2．37a）式の右辺最終項に現れる衝突項の非線形項 Cab（ fˉa1，fˉb1）についても考慮したドリフト

が磁気面関数ではなくなる．この偏りの影響を含めた計
算は炉心プラズマに入ったタングステンの輸送を研究する

運動論方程式ソルバーも開発されており［34］，これは新上で重要と考えられ，近年盛んに研究されている
［37，
38］
．

古典輸送シミュレーションとしては最も複雑かつ大規模さて，（3．1．2）式では平衡との観点から速いオーダリ
ングのフローは軸対称系においてのみ考えられるとしたなコードの部類に入る．
ここまで，線形化衝突項の扱いについていくつかのタが，ここで考えている速い平均流 u0 は径電場による E×B

イプを見てきたが，他にもまだ衝突項の拡張が必要とな

フロー由来であり，粒子種に依らず共通である．一方，る場合がある．それはプラズマが外部から NBI などに

熱速度 vth,a∝ Ta/ma は粒子種に依存し，熱核融合炉ではよって強く加熱されている場合や，将来の核融合炉の DT

イオン ‐ イオン間衝突による温度緩和によりどのイオン反応によって生じる高速 α 粒子によるバルクプラズマの自種もほぼ等温であると考えられるから，トカマク以外の

己加熱を扱う場合である．このような場合，加熱を引き

非軸対称系でも，特に両極性条件を満たす径電場が大き起こす高速粒子の分布は Maxwell 分布をゼロ次としてそい場合には，バルクの水素や電子に対しては u0/vth ＜ δ1

こからのずれを 1 次とする，という線形化衝突項での扱であっても，タングステンのような重イオンに対しては

いにそぐわない状態になっている．こうした場合にどの

u0,I/vth,I ～ δ 0 という状態に容易になると想定できる．したように衝突項を扱うかについては，3．2 節で具体的な事例

がって，速いフローのオーダリングを用いた新古典輸送を紹介したい．

計算法は，ヘリカル型核融合炉の不純物輸送問題を扱う

5 ）速い平均流の存在下でのドリフト運動論方程式：本上でも重要である．しかしながら，速いフローのオーダ

節ではこれまで，（3．1．1）式をベースに新古典輸送計算にリングのドリフト運動論方程式はこれまで軸対称系に対おける様々なレベルの近似について述べてきた．今後核してのみ定式化とコード開発が進められてきており，非融合炉の運転予測に新古典輸送計算を利用するにあたっ

軸対称系への拡張は今後の課題である．

て，もう 1 つドリフト運動論方程式の拡張が必要になる

本節の最後に，新古典輸送計算法の計算モデルの差に点がある．ここまでのドリフト運動論方程式では，ドリ

よってどのように計算結果が変わるか，いくつか事例を

フトオーダリングにおける暗黙の仮定としてプラズマの紹介する．一つ目は，ヘリカルプラズマにおける両極性

平均流は uˉ a ～ O（δvth,a）程度と仮定されていた．しかし例

径電場の計算における径方向局所近似と有限軌道幅効

えば軸対称トカマクにおいては，トロイダル回転を減衰

果を含めた大域的計算の比較である．図 1 は，W7-X のさせる機構は O（δ1）次までのモーメンタムバランスの（2．

ECH 加熱放電において，観測された温度密度分布から新

古典輸送計算を行い，両極性条件 ΓiNC＝ΓeNC を満たす径電21）式には存在しないため，接線 NBI 等外部からのトルク
入射があると容易にトロイダル方向に回転する．今，準場 Er の予測をいくつかのシミュレーションで行い比較し対称系配位まで話を拡張し，eζ－方向に磁気座標系で見たたものである［39］．この中で，DKES と SFINCS は E×B

磁場の対称性があると仮定しよう．今まで O（δ0）次と仮定

ドリフトのみを入れた径方向局所近似，FORTEC-3D はしていたポテンシャル Φ0 が eaΦ0/T0～ O（δ－1）の大きさを

磁気ドリフトの磁気面接線，垂直の両成分を入れた有限持つとする．この場合，平均流として E0 ＋ u0 × B ＝ 0 を満たす E×B トロイダルフロー u0 ＝ uζ0eζ，uζ0 ＝－ Φ0' /ι が存在する．この平均流を含めた MHD 平衡の force balanceの式は以下のようになる［35］．

（∑ m n ）u ・∇u ＝J×B－∇p ,1∂g＝ιJ ． （3．1．2）（上式・e ）→ （∑ m n （u
） ）  
∂ζ
2
a

a a

ζ

0

0

a

0

a a

ζ 2
0

ζζ

s

ここで gζζ ＝ eζ・eζ．軸対称トカマクの場合，∂gζζ/∂ζ ＝ 0
である．一方，準対称配位では一般的に ∂gζζ/∂ζ≠0 であるから，（3．1．2）式から MHD 平衡を支えるのに有限の径

図1

方向電流 Js が必要となる．しかし，磁気座標系を構築す
317

W7-X の ECH 加熱放電における，両極性径電場の予測計
算の計算モデル間の比較［39］.

Journal of Plasma and Fusion Research

Vol.100, No.7

July 2024

軌道幅効果を含む大域的計算となっている．このプラズの新古典理論ではつなぎ変えの起こる磁気面を決定でき

マでは，実験観測でも磁気軸側が正，周辺部で負の径電ないことと矛盾しないが，大域的新古典シミュレーショ

場が観測されているが， 3 つのシミュレーションともほぼ

ンで Γanom を入れなくても滑らかに正負の電場をつなぐ解

同じ傾向を示しており，この図には示されていないものが得られることとは一致しない．また，Maxwell 構成則の，値としても観測値に近いものとなっている．局所近では電場の遷移はある磁気面で不連続に起こることにな

似と大域的計算の差は，正負の電場の切り替わる磁気面

るが，実験観測から径電場はある径方向の幅で連続的に近くに現れている．2．4 節と本節の 1 ）で説明したように，
FORTEC-3D ではプラズマ全域で径電場の時間発展込み

変化していると見られている．局所近似の新古典フラックスに，有限軌道幅効果の一端として電場のシアー ∝Φ''でドリフト運動論方程式を解いている．その結果，正負のに比例する項をモデル的に入れることで［41］ある程度径

径電場の入れ替わる付近でも連続的な径電場の遷移が自方向に幅のある遷移領域を持った両極性径電場分布の予然に求められる．一方，局所近似計算では磁気面毎に独測がされているが［42］，径電場シアー効果のモデルは多立した計算となっており，径電場をパラメータとして振分に定性的であり，定量的な正確さには疑問が持たれるりながら両極性条件を満たす Er を探索する．この図でもが，理論的進展がほとんどされていない．大域的シミュ

一部，正負 2 通りの解が存在している領域があるが，そレーションでの径電場の遷移面付近では，この径電場シのどちらが実際の解として選ばれるか，どこでどのようアーの効果の他，解析的モデルでは入っていない径方向に電場が遷移するかは，局所近似計算だけで決定するこのモーメンタム輸送などの効果も影響していると考えらとはできない．新古典輸送と乱流輸送の相関を考える場れる．また， 2 ）で取り上げた磁気面上のポテンシャル変合，例えば径電場シアーによる乱流の抑制効果を考える

位 Φ1 は，径電場の正負で位相が逆転するため，これも間

上では，この径電場のつなぎ変わり面でどの程度の幅で接的に径電場遷移の位置決定に関係している可能性があ

正負の径電場が遷移しているかが重要となると考えられる．径電場分岐の問題は今後のより一層の研究の発展がる．一方，それ以外の位置では 3 者の計算結果がよく一致待たれる．

している．これは，W7-X が新古典輸送を抑えるために磁2 つ目は， 3 ）で取り上げたいくつかの径方向局所近似

気ドリフトによるドリフト軌道の磁気面からのずれを抑の比較研究である．［12］では，局所近似版の FORTEC-

制するような磁場配位を持つように設計されているため，

3D においていくつかの異なる径方向局所近似の計算を

有限軌道幅効果があまり径方向輸送に影響しないためで比較している． 1 つ目は径電場による E×B ドリフトのあると考えられる．このような観点から，準対称配位なみ残し，磁気ドリフト項を完全に落としたモデル（Zero-

どのいわゆる新古典最適化配位に近い場合，径方向局所

Magnetic-Drift, ZMD）， 2 つ目は ZMD から更に vE×B ＝モデルは定量的にも十分な精度の計算が可能であると言E×B/B2 → E×B/
〈B2〉と E×B ド リ フ ト の 非 圧 縮 近 似 をえる．なお，大域的計算法による時間発展計算においては，行ったもの（DKES-like）， 3 つ目は ZMD に磁気ドリフ

複数の両極性解が存在している場合も一方の解に自然にトの磁気面接線成分を足したモデル（Zero-Orbit-Width,

径電場が収束していくが，この場合もなぜ一方の解が選ZOW）である．図 2 に LHD 磁場配位で径電場を振った

ばれたのかについては説明することがまだできていない．

場合のそれらのモデルでの新古典粒子イオンフラックス

従来の理論では，ある磁気面 s において，正負 2 つの両極の依存性の比較を行っている．また，バウンス平均計算性条件を満たす径電場の根 E（s）
＝ E＋（s），E－（s）があるr

コードの GSRAKE と，近似のない大域的なドリフト軌道

場合に，Maxwell 構成則（Maxwell’s construction）といを解く Global 版 FORTEC-3D の計算結果とも比較していう考え方に従って，る．径電場ゼロ近くで磁気ドリフト項が入っていないモ

  E＋



ΔΦ（s）≡

E－

NC
NC
dE（Γ
＝0，
r
i －Γe ）

デル（ZMD, DKES-like, GSRAKE）の Γi が発散傾向にあるのに比較し，ZOW ではピークが負電場側にシフトし，（3．1．3）

またそのピークも大幅に下がっている．ZOW の計算は有すなわち電場が二つの根の間でつなぎ変わっても径方向

限軌道幅効果の入った global 計算の結果を完全には再現に電流が流れないところで遷移が起こる，というモデル

しないもののその傾向をよく捉えており，磁気ドリフト

NCが提案されている［40］．ただしこのモデルでは，Γi,eと項の接線成分が径電場ゼロ付近でのヘリカル系の新古典

して局所近似の解析解を使い，また（2．34）式のように駆輸送を正しく評価する上で重要となることを示した結果

pa'
動 力 Xa1＝－ －eaΦ'，Xa2＝－Ta' を 与 え て，Xa1 の 中 のnaとなっている．また，図 2 では径電場が大きい領域での比較も示されており，E×B ドリフトの非圧縮近似を使っΦ'∝－Er を変化させるのではなく，乱流輸送など新古典ている DKES-like モデルでは，理論的に予測されている輸送以外のプロセス（ただし両極性な）を含めたトータ
ルフラックス Γt ＝ ΓNC ＋ Γanom，Γanom∝n' を想定し，（3．

Mp≃1 付近に現れる Γi のピークが再現できないことが確

1．
3）式の積分中で Er を変化させる際に Γt を保存するように n' を変化させるという考え方をする．このモデルでは，軌道に対する近似は新古典輸送の評価に大きく影響を与認できる．このように，局所近似計算におけるドリフト
えることが発見された．なおこのような磁気ドリフト項

Γanom＝0 としてしまうとそもそも電場のつなぎ変えが起の接線成分の効果は SFINCS，KNOSOS コードでも確認されている．なお［12］では電子フラックスについてもこらない，という結果になってしまう．これは局所近似
318

Special Topic Article

3.1 Numerical Schemes for Drift-Kinetic Equation and Several Approximations in them	S. Satake

参考文献

図2

［ 1 ］J.M. Garca-Regaa et al., Nucl. Fusion 57, 056004(2017).
［ 2 ］K. Fujita et al., J. Plasma Phys. 86, 905860319(2020).
［ 3 ］S. Satake et al., Plasma Fusion Res 1, 002 (2006).
［ 4 ］R. Kleiber et al., Comput. Phys. Commun. 182, 1005
(2011).
［ 5 ］A. Molln et al., Plasma Phys. Control. Fusion 60,
084001 (2018).
［ 6 ］H. Yamaguchi and S. Murakami, Nucl. Fusion 58,
016029 (2018).
［ 7 ］S.P. Hirshman et al., Phys. Fluids 29, 2951 (1986).
［ 8 ］W.I. van Rij and S. P. Hirshman, Phys. Fluids B 1,
563 (1989)［ 9 ］D.A. Spong, Phys. Plasmas 12, 056114 (2005).
［10］C.D. Beidler and H. Maaßberg, Plasma Phys. Control. Fusion 43, 1131 (2001).
［11］S. Matsuoka et al., Phys. Plasmas 22, 072511 (2015).
［12］B. Huang et al., Phys. Plasmas 24, 022503 (2017).
［13］J.L. Velasco et al., J. Comput. Phys. 418, 109512(2020).
［14］M. Honda et al., Nucl. Fusion 52, 023021 (2012).
［15］H. Sugama et al., Phys. Plasmas 23, 042502 (2016).
［16］M. Landreman et al., Plasma Phys. Control. Fusion
56, 045005 (2014).
［17］V. d’Herbemont et al., J. Plasma Phys. 88, 905880507(2022).
［18］S. Kitajima et al., Nucl. Fusion 53, 073014 (2013).
［19］K.C. Shaing, Phys. Fluids B 5, 3841 (1993).
［20］N.C. Logan et al., Phys. Plasmas 20, 122507 (2013).
［21］K.C. Shaing et al., Nucl. Fusion 50, 025022 (2010).
［22］H. Sugama and S. Nishimura, Phys. Plasmas 9, 4637
(2002).
［23］H. Sugama and S. Nishimura, Phys. Plasmas 15,042502 (2008).
［24］A.H. Boozer and G. K. Petravic, Phys. Fluids 24, 851(1981).
［25］M.N. Rosenbluth et al., Phys. Rev. 107, 1 (1957).
［26］S. Satake et al., Comput. Phys. Commun. 181, 1069
(2010).
［27］R. Kleiber et al., Comput. Phys. Commun 295,
109013 (2024).
［28］H. Sugama et al., Phys. Plasmas 16, 112503 (2009).
［29］S. Satake et al., Comput. Phys. Commun 255, 107249(2020).
［30］E.A. Belli and J. Candy, Plasma Phys. Control. Fusion 54, 015015 (2012).
［31］M. Landreman and D. R. Ernst, J. Comput. Phys.
243, 130 (2013).
［32］H. Sugama et al., Phys. Plasmas 26, 102108 (2019).
［33］S. Matsuoka et al., Phys. Plasmas 28, 064501 (2021).
［34］R. Hager et al., J. Comput. Phys. 315, 644 (2016).
［35］H. Sugama et al., Plasma Phys. Control. Fusion 53,
024004 (2011).
［36］F.L. Hinton and S. K. Wong, Phys. Fluids 28, 3082
(1985).
［37］C. Angioni and P. Helander, Plasma Phys. Control.
Fusion 56, 124001 (2014).

 々な径方向局所近似モデルと大域的モデルによる，LHD
様
磁場配位での新古典イオンフラックスの径電場依存性の比較［12］．上図が径電場ゼロ付近，下図が径電場の大きい領域の比較．

同様の比較を行っているが，磁気ドリフト項を無視した
影響はイオンにのみ顕著に表れ，電子に対しては ZMD，
DKES-like モデルでも大域的計算と遜色ない計算ができることがわかった．径電場ゼロ付近で磁気ドリフト項を無視すると Γi が過大に評価されてしまうことは，いわゆ
る 1/ν 領域の新古典輸送の解析解に相当するものをシミュ
レーションで見ていることに相当する．この径電場ゼロ付近の偽のピークによって，場合によっては存在しない
両極性条件 Γi ＝ Γe を拾ってしまうことになるので注意が必要である．また［12］では，衝突項のモーメンタム保存性の有無による新古典輸送計算の結果の違いを準ヘリカ
ル対称配位 HSX で比較もしているので，興味のある方は読んでみてほしい．
以上の例に示したように，新古典輸送計算は採用する近似によって計算結果が異なる場合があるので利用するにあたってそのコードの特性を理解することが重要である．本稿が核融合研究者の間でそのような意識づけのきっかけになれば幸いである．

謝

辞

本章にまとめた内容は筆者が総合研究大学院大学で博士課程の指導教員を務めた松岡清吉氏，Huang Botsz 氏，藤田慶二氏の学位論文における研究成果によるところが大きい．この場をお借りして改めて御三名の新古典輸送
理論研究への貢献に感謝を申し上げます．本章の執筆は科学研究費［基盤研究（C）21K03517］の援助を得て行わ
れました．
319

Journal of Plasma and Fusion Research［38］D. Fajardo et al., Plasma Phys. Control. Fusion 65,
035021 (2023).
［39］N.A. Pablant et al., Phys. Plasmas 25, 022508 (2018).
［40］K. Itoh et al., J. Phys. Soc. Jpn. 70, 1575 (2001).

Vol.100, No.7

July 2024［41］D.E. Hastings, Phys. Fluids 28, 334 (1985).
［42］S. Toda and K. Itoh, Plasma Phys. Control. Fusion
53, 11501 (2011).

320

J. Plasma Fusion Res. Vol.100, No.7 (2024) 321-327

小特集 古くて新しい新古典輸送理論の新展開

3 ．新古典輸送計算の実際～近似モデルと新展開～

3. Neoclassical Transport Calculations in Practice:
Approximate Models and New Developments

3．2

高エネルギー粒子を含む問題における Coulomb 衝突項

3.2 The Coulomb Collision Operator for Problems Including Energetic Particles
西村伸，奴 賀 秀 男

NISHIMURA Shin and NUGA Hideo
核融合科学研究所

構造形成・持続性ユニット（原稿受付：2024年 4 月25日）

加熱されたプラズマにおける熱化粒子挙動を考える時のために，Maxwell 分布から著しくかけ離れた高エネ
ルギー粒子速度分布も場の粒子種として扱う Rosenbluth ポテンシャルの取り扱い方法を示す．また非線形衝突効果の一例として，過去の高エネルギー粒子関係理論でしばしば無視されがちであった高エネルギー粒子どうしの衝突による速度分布変化を示唆する中性粒子ビーム入射（NBI）加熱実験を紹介する．
Keywords:
Coulomb collision, NBI-heated and/or burning plasmas, energetic particles, Rosenbluth potentials, nonlinear collision
operator

3．2．1

はじめに

しかし今後の ITER では中性粒子ビーム入射（NBI）によ

衝突項の取り扱いについて様々な工夫がなされてきたる高速イオンにより核融合生成高速イオンの速度分布が歴史経緯は前章 3．
1 に述べられた通りであるが，これらは非等方になる可能性が指摘されている．そのような高速

基本的に熱化粒子どうしの衝突の話である．今回担当編

イオンどうしの衝突による速度分布変化を示唆する実験

集委員と企画取りまとめ責任者より，加熱や核融合反応

結果が複数 NB を用いた既存装置実験でも得られつつあるで生じた高エネルギー粒子を含む状態における衝突項に事を紹介する．

関する解説を依頼された．高エネルギー粒子が関わる問

3．2．2

題と言っても，テスト粒子をそれら高エネルギー粒子とし，場の粒子として Maxwell 分布となっている熱化粒子

 熱されたプラズマにおける熱化粒子衝
加突項を想定する線形化衝突項に基づいて高エネルギー粒子速これは「古くて新しい」というこの小特集タイトル

度分布を決定する理論やシミュレーション技法についての「古い」話である．NBI 加熱プラズマや燃焼プラズマは過去より多数の論文が出版されている．この章ではそのような高速イオンが存在する状態における熱化粒子のれらであまり論じられなかった二つの問題を取り上げる．

挙動を考えるには，高速イオン減速速度分布関数のよう一つはそのような高エネルギー粒子の存在下での熱化粒に Maxwell 分布から著しくかけ離れた速度分布を場の粒子挙動を考える熱化粒子運動論的方程式側における衝突子として取り扱う Rosenbluth ポテンシャルの計算方法が項である．もう一つは高速エネルギー粒子どうしの衝突必要である．そのような計算式はジャイロ位相平均速度である非線形衝突効果である．テスト粒子種と場の粒子分布については論文［1］に Rosenbluth, MacDonald, Judd

種の速度分布両方が著しく Maxwell 分布からかけ離れて（以下では RMJ という）が既に示していたのである．しいる衝突項は非線形衝突項とも呼ばれ，3．1 章で説明さか し こ の 式 は，Review of Plasma Physics な ど を 始 めれたように熱化粒子どうしの衝突ではそのような非線形とするその後の Coulomb 衝突項に関する解説書籍や新項を落とす線形化近似を用いるのが通常である．高エネ

古典理論レビュー論文などに取り上げられる事がほとん

ルギー粒子関係理論でも，そのような高エネルギー成分

ど無く，非線形衝突項を考える研究者や，熱化粒子同士f（x，v）が粒子数としては熱化粒子より少ない状態を想定
fを扱う線形化衝突項でも場の粒子の速度分布として任意
Legendre 展開項数を取り入れようとする研究者（論文［2］ff）を落とすことがしばしばであった．
して非線形項 C（fff f，National Institute for Fusion Science, Toki, GIFU 509-5292, Japancorresponding author’s e-mail: nishimura.shin@nifs.ac.jp

321Ⓒ2024 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research

Vol.100, No.7

July 2024など）にしか知られていなかった．こうなったのは古典

ム駆動電流［4］，燃焼炉心のブートストラップ電流［5］な及び新古典輸送理論の歴史が関係しており，場の粒子項

ど，高速イオン存在下でのプラズマフローを考える時にCab（faM，fb）を含む問題の取り扱い手法として Braginskii

は場の粒子として Maxwell 分布から著しくかけ離れた高の行列表現［3］が普及したからである．Braginskii がこの

速イオン減速速度分布を扱うことになるので本来は使え行列要素を求める方法として示したのは Landau 型表現ない．しかしここでも論文［1］のようなポテンシャル計算に Sonine 多 項 式（ 一 般 化 Laguerre 多 項 式 ）L（α）
≡
j （K）
K －α
j －K j＋αj式が必要であると広く認識されるには至らなかった．そ）/dK の母関数を代入する方法であ

れは高速イオン駆動の電流を目的とする場合，熱化イオ

り，Rosenbluth ポテンシャルという概念を経ずに行列要

ンフローを無視して電子フローだけ計算することでも一L（l＋1/2）
L（l＋1/2）
素 vlP（ξ）（xa2）Cab（vlP（ξ）（xa2）faM，fbM）d3v，ll
jk

応用が足りて，その目的に電子 - イオン衝突に関する小質（e K/j!）d（e K＋ 1/2） 2
L（lL（l＋1/2）
v P（ξ）（x a ）C ab（f aM ，v l P（ξ）
（x b2 ）f bM ）d 3 v
l
l
j
k
l
l
l
－1 l
2
［P（ξ）
≡（2 l!） d（ξ －1）/dξ：Legendre 多項式，ξ ≡ v‖/v，l2
xa2 ≡ mav2/2Ta ≡ v2/vTa（－xa2）：，faM ≡ naπ － 3/2v－3
Ta exp量比（small-mass-ratio）近似が使えてしまったからであlる．1970〜1990 年台に高速イオンに関わる解析的理論を示した諸論文に vb2≪2Te/me の仮定がしばしば登場する一Maxwell 速度分布］が得られる．これを用いて運動論的つの理由はこれである．しかしその後荷電交換分光がビー問題を代数方程式に変換しようという考えは「新古典」

ム駆動のイオンフローを調べる目的に使われるようにな

以前の「古典」の時代に現れたが，「新古典」の時代にると［6］これだけでは済まないこととなった．更に将来のなってからの 1970 年台に Hirshman と Sigmar が複数イ

準対称配位実験ではこのフロー発生のみならずそれに伴オン種状態の取り扱いのためにこれを積極的に応用し始うビーム駆動の径方向輸送［7，8］も調べられることになろめ，その集大成が有名な論文［4］である．この代数的取りう．ここでは論文［1］の式をジャイロ位相平均に限らない扱いはあくまでもテスト粒子種と場の粒子種が Maxwell

一般的速度分布に拡張する形で球座標計算式を示すとと分布に近づきそれぞれのエネルギー空間構造が有限項

もに，その幾つかの応用を示す．なお，本 3．2 節は CGS-

数 Laguerre 展開で近似できる場合の手法である．ビー

Gauss 系で記述する．

3．2．2．1 一般的速度分布に対するRosenbluthポテンシャル
次の Rosenbluth ポテンシャル定義式から出発する．
f（v'）b
H（fb）≡     d3v'，G（ fb）≡ f（v'）|v－v'|d3v'
b|v－v'|



（3．2．1）

速度分布関数は球座標形式の速度空間座標（v，θ，ϕ）における次の球面調和関数による展開形で与えられると想定する．
∞

l

l＝0

m＝1

0
P（cos
cos（mϕ）＋bm
sin（mϕ）}
f（v）＝  
θ）＋ ∑ Plm（cos θ）{ambl∑ a（v）
l
l（v）
l（v）
∞

θ，ϕ）
≡∑ f（l）
b （v，（3．2．2）l＝0dξm：Legendre 陪関数］次に v，v' の間のコサインを z とした時の［Plm（ξ）≡（1 － ξ2）m/2 dmP（ξ）
l|v－v'|±1＝（v2－2vv'z＋v ' 2）±1/2，z＝cos θ cos θ '＋sin θ sin θ ' cos（ϕ－ϕ'）ν（ν ＝ 1/2 の 時 に Cn1/2（z）＝ P（z）） の 母 関 数に は 文 献［9］の ア イ デ ア を 応 用 し，Gegenbauer 関 数 C（z）
n
n
∞

ν
－ν
hn（|h| ＜ 1）を使って Legendre 多項式 P（z）
（1－2hz＋h2）＝  Cによる展開表現とする．|v － v'|－1 はこの母関数
n
∑ （z）
n
n＝0

その物であって，│－1|v－v'| ＝
│


（ ）
v
v'
1
  P（z）（ v' ） for v ＞ 1
v' ∑
v' n
v'
1 ∞
for ＜ 1  P（z）
n
∑
vn＝0
v
v
n＝0

（3．2．3）n∞nである．次に |v － v'| はやや複雑な導出手順をとるが，その結論を示すと

（ ）
（ ）

（ ）
（ ）

n
∞
v' 2
v'
1 ⎫
1
 v   P（z） v' ⎧⎨   －   ⎬ for
＜1
n
∑
│ n＝0
v ⎩ 2n＋3 v
v
2n－1 ⎭|v－v'|＝
∞v n⎧ 1
v 2
v'
1 ⎫
│
⎨   
  P（z）－    ⎬ for＞1n∑ v' n＝0v' ⎩ 2n＋3 v'
v
2n－1 ⎭

322（3．2．4）

Special Topic Article

3.2 The Coulomb Collision Operator for Problems Including Energetic Particles	

S. Nishimura and H. Nugaとなる．ここでn!
P（z）＝P（cosθ）P（cosθ '）＋2 ∑（n－m）    Pnm（cos θ）Pnm（cos θ '）{cos（mϕ）cos（mϕ'）＋sin（mϕ）sin（mϕ '）}nnnm＝1（n＋m）!

（3．2．5）

であること［10］に注意する．式（3．2．2－5）を式（3．2．1）に代入し，球面調和関数の直交関係を利用した d3v' 積分を遂行すると以下の諸式が得られる．
  

  

∞
v
∞ f（l）1 ⎧ 1θ，ϕ） ⎫b （v'，
l＋2 （l）
H
fb （v'，θ，（fb）＝4π ∑    ⎨  ϕ）dv'＋vl      dv' ⎬
l＋1 （v'）
l－1
⎭
（v'）
v
0
l＝0 2l＋1 ⎩ v

（3．2．6）

  

  

∞
v
∞ f（l）∂H（ fb）
1 ⎧θ，ϕ） ⎫1b （v'，
l＋2 （l）
dv' 
fb （v'，θ，ϕ）dv'＋lvl－1     ⎬    ＝4π ∑     ⎨－  
l－1
l＋2 （v'）
⎭
v
∂v
（v'）
v
0
l＝0 2l＋1 ⎩

（3．2．7）

  

  

∞
v
l－1 1 v∂G（ fb）
1 ⎧ l＋1 1
l＋4 （l）
fb （v'，θ，ϕ）dv'＋    l （v'）l＋2 f（l）
θ，ϕ）dv'    ＝4π ∑      ⎨－     b （v'，
l＋2 （v'）∂v
2l＋3 v
2l－1 v 00
l＝0 2l＋1 ⎩  ∞ （l）  ∞ （l）
l＋2
fb （v'，θ，ϕ）
fb （v'，θ，ϕ） ⎫1dv'－   vl－1
dv' ⎬
     
     
            ＋   vl＋1
l－1
l－3
⎭
（v'）
2l－1
（v'）
2l＋3
v
v








  

（3．2．8）  ∞vv（l＋2） 1（l－1）l 1∂2G（ fb）
1 ⎧（l＋1）
l＋4 （l）
fb （v'，θ，ϕ）dv'－      （v'）l＋2 f（l）
θ，ϕ）dv'    ＝4π      ⎨       
  ∑b （v'，
l＋3 （v'）l＋12
v
2l－1 v2l＋3∂v
00
l＝0 2l＋1 ⎩（3．2．9）  ∞ （l）  
fb （v'，θ，ϕ）
（l＋1）
（l＋2）
（l－1）l l－2 ∞ f（l）
θ，ϕ） ⎫b （v'，
dv'－    v
dv'
     
     
            ＋       vl
⎬
l－1
l－3
⎭
2l－1
2l＋3
（v'）
（v'）
v
v








これは論文［1］の積分公式を任意の球面調和関数（l，m）を含められるように拡張した物となっている．
3．2．2．2

応用例 1: 線形化衝突項の場の粒子項の計算

この一つの応用先は線形化衝突項の内の，場の粒子項

（ ）

（

）

eaeb 2ma ⎧4πTama∂H（ fb） mav2 ∂2G（ fb）⎫⎨   fb－H（ fb）＋
Cab（ faM，fb）＝4π    ln Λab faM －1 v    ＋      ⎬maTa ⎩ m b
mb
2Ta
∂v
∂v2 ⎭

（3．2．10）（ fb），
である．この表式は RMJ 表現に Maxwell 分布 fa ＝ faM を代入した上で Poisson 方程式 ∇v2 H（ fb）＝－4π fb，∇v2 G（ fb）＝2Hv）の存在下で Caf（ faM，ff）を含めた熱化粒子のドリフト運動論方∇v2 ≡∑ ∂2/∂vα2 を用いて得られる．高速イオン速度分布 f（x，fα

程式（DKE）の解を求めるには，上記に述べたように当初 Hirshman-Sigmar に提唱されその後一般的非対称トーラスにも拡張された，ピッチ角空間構造 Legendre 展開成分（特にその一次）のエネルギー空間構造に関する Laguerre 展開係数の従う連立代数方程式を導く手順（本章論文［7］及び 3．1 章の論文［9，22-23］）を取る．このために DKE 項の Laguerre 展開係数が必要である．Cab（ faM，fb）にそれを行う時は，不定積分公式
exp（－a2x2） nn!

      x∑（n－k）
x exp（－a x ）dx＝－      !a2ax（2n－1）!!
（2n－1）!! πexp           ＋       Φ（ax）（－a x ）∑
x exp（－a x ）dx＝－    
2 a （2n－2k－1）!!
2a
2 a
2n＋1

2 2

2n

2 2

2

2

k＝0

2 2

2k

2n－2k－1n－1

k＝0（3．2．11）2（n－k）k 2kn＋1 2n＋1n－2
x2n－2k－1（2n－1）!!
（2n－1）!! π 2
exp（－a2x 2）∑
G（ax）．
          ＋       x＝－    k 2k2!!
2a2na2n－1
k＝0 2 a （2n－2k－1）（3．2．12）

を用いて dv の部分積分を行う．ここで
  

∞
x
2
2
（－1）n
2
       x 2n＋1
Φ（x）≡   e－y dy＝   ∑π 0
π n＝0（2n＋1）n!

∞
Φ（x）－x Φ'（x） 2（－1）nG（x）≡       ＝         x 2n＋1∑22xπ n＝0（2n＋3）n!

はそれぞれ誤差関数と Chandrasekhar 関数である．紙面の都合であまり広範囲の Legendre 展開次数，Laguerre 展開次数について Cab（ faM，fb）展開係数を示す事はできないが，三例だけ示しておこう．より詳しくは論文［11］を参照されたい．



ma v2Cab（ faM，fb）d3v＝－mb v2Cba（ fb，faM）d3v（ ）

（ f dξ）v dv

  

1/2 ∞
2
n（e
ma 1
⎧
⎫
a aeb）ln Λab ma⎨xaG（xa）－  exp（－xa2）⎬＝32π2        
⎭
ma
mb π
2Ta
0 ⎩

323

  1

－1

b

2

（3．2．13）

Journal of Plasma and Fusion Research

Vol.100, No.7

July 2024（ ） （ ξf dξ）v dv （3．2．14）n（e e ）ln Λx⎧ m⎫
f dξ v dv （3．2．15）m vξL （x ）C （ f ，f ）d v＝24π         ⎨（－x ）－G（x ）（⎬  ξ
（ m ＋1）  exp
）
⎩
⎭
T
π
2
n（e
a aeb）ln Λab ma
ma vξCab（ faM，fb）d3v＝－mb vξCba（ fb，faM）d3v＝8π2        ＋1Ta
mb（3/2） 2
a
1

a

fb ≡
［

ab

aM

3

b

a

2

2

a b

ab

  ∞

a

0

a

  ∞

  1

0

－1

a

G（xa）

  1

2
a

b

2

b

a

－1

2

b

  π

 f dϕ/（2π）］ 式（3．2．13－14）においてはエネルギー・運動量保存のために場の粒子項 C （ f ，f ）とテスト粒子項
－π

b

ab

aM

b

Cba（ fb，faM）の積分が一致しているが，テスト粒子項のこの種の積分手順については論文［12］を参照されたい．このようなL（l＋1/2）
積分公式に fb＝vlP（ξ）（xb2 ）fbM を代入すると論文［3，4］などに示されてきた Braginskii 行列要素も再現する．これにljは部分積分で得られる漸化式
  ∞

x
0

と

  ∞

（

  ∞

0

1

 x Φ（x ）exp（－x ）d x ＝ 2 {1＋（v /v ）}
a

0

2
a

b

a

Tb

2 －1/2Ta）（2n＋1）/2
2 －
vTa

1＋ 2
vTb

vTa（2n－1）!!
2n－1
Φ（xb）exp（－xa2 ）dxa＋
    
a
vTb
2n＋1

x

2n＋1
Φ（xb）exp（－xa2 ）d xa＝na（3．2．16）を使う．また式（3．2．6－12）は，場の粒子種 b が熱化粒子である場合にはL（l＋1/2）
Cab（ faM，vlP（ξ）（xb2 ）fbM）を知る目的にも使え，これも Braginskii［3］が示した行列要素導出手順では得られない事で
ljある．

3．2．2．3

応用例 2 ：非線形衝突項による粒子種間エネルギー交換

もう一つの応用は非線形衝突項による粒子種間エネルギー交換である．非線形衝突項それ自体は複雑だが（線形化衝突
項は Landau もしくは RMJ のデカルト座標表現を球座標にする速度空間変数の変換で取り扱いやすくなるが，非線形のままでは座標の工夫をしても簡単にならない），そのエネルギー交換積分は簡単な式に帰着する．デカルト座標での部分積分を行い，Poisson 方程式 ∇v2 G（ fb）＝ 2H（ fb）を使うと（ ）  （）
eem∂H（ f ）⎫⎧ln Λ f ⎨H（ f ）＋ 1＋
⎬d v
＝8π
（  
（ m ）v   
⎩
m ）
∂v ⎭
ee
1
ln Λ ∑
     
＝2
（4π  
m ）
2l＋1
f （v'，θ，ϕ） ⎫mm1⎧
l＋l＋1 v        dv'⎬d v
  ×f ⎨－（ （l＋1）＋l）  （v'） f （v'，θ，ϕ）dv'＋（） （v'）⎩⎭vmm
eaeb 2∂ ma ∂H（ fb） 1 ∂2G（ fb） ∂fa 3
dv
d3v＝－4π    ln Λab v2  
v2Cab（ fa，fb）
fa    －     ma∂vα mb∂vα
2 ∂vα∂vβ ∂vβ



a b

2

ab

a

a b

2

a

a

b

b

3

b

∞

ab

a

l＝0

  v

a

a

l＋1

b

a

l＋2 （l）b0が得られる．球面調和関数の直交性を使いl  ∞ （l）
b
v

b

3

l－1

  ∞

 dv に関する部分積分を行うとこの式は0ln Λ ∑     
v C （ f ，f ）d v＝2（4π   
m ）
2l＋1
m
f （v'，θ，ϕ） ⎫1⎧ m⎨l＋l＋1   （v'） f （v'，θ，
ϕ）dv'－（ （l＋1）＋l）v        dv'⎬f d v
×（
）v
⎩ m
⎭
m
（v'）
2

ab

a

a
b

b

eaeb

3

2

  v

l＋1

0

∞

ab

a

（3．2．17）

1

l＝0

a

l＋2 （l）al  ∞ （l）
a

v

b

（l） 3
b

l－1

（3．2．18）2（ fa，fb）＋ mbCba（fb，fa）］d3v ＝ 0 に対応している．fa ＝ faM もしくとも書き換えられる．この事はエネルギー保存 v［maCab

は fb ＝ fbM を代入すれば線形化衝突項のエネルギー交換式（3．2．
13）になることは言うまでもない．この fa，fb の片方が熱化
粒子速度分布であって Legendre-Laguerre 展開形で与えられているなら，それを含む

  v

  ∞

0

v

 dv'， dv' の不定積分には式（3．2．

11-12）を用いる．更にもう片方も熱化粒子であるなら式（3．2．16）を用いて Braginskii 行列要素と同様に熱速度比 vTa/vTb の関数としての定積分値を得る．

3．2．3

非線形 Coulomb 衝突項

3．2．3．1

Coulomb 衝突の非線形性

ネルギー粒子の速度分布関数を評価する場合において，Coulomb 衝突の非線形性が考慮される．加熱以外では，

Coulomb 衝 突 の 非 線 形 性 と は， 速 度 分 布 関 数 の 非トカマクにおけるディスラプション時の逃走電子生成過Maxwell 成分同士の Coulomb 衝突のことであり，速度分程が例として挙げられる［13］．

布関数の非 Maxwell 成分が大きい場合に顕著に現れる．

例えば，中性粒子ビーム（NBI）によって生成される

このような状況は，外部加熱などによって，大量の高エNB 高速イオンの速度分布関数について考える．このとネルギー粒子が供給される場合に生じるため，主に高エき，ビームのガス種と，プラズマのイオン種は同じもので
324

Special Topic Article

3.2 The Coulomb Collision Operator for Problems Including Energetic Particles	

S. Nishimura and H. Nuga

あるとする．NB 高速イオンは背景プラズマとの Coulomb

による損失等を考慮しない場合）．しかしながら，3.5 MeV

衝突によって，減速・ピッチ角散乱される．通常，このの α 粒子と 1 MeV の重水素 NBI の速度比は vHe/vD ～1.3 程過程は背景プラズマを構成する電子・イオンの速度分布度であるため，ビーム方向のピッチ角においては相対速を Maxwell 分布であると仮定して Fokker-Planck 方程式

度が非常に小さくなる．このため，ビーム方向のピッチを用いて計算される．しかし実際には，NB 高速イオン自角において，重水素ビームイオンと高エネルギー α 粒子と身によって，イオンの速度分布関数は非 Maxwell な分布の Coulomb 衝突が無視できなくなり，結果として，α 粒へと変形している．NB 高速イオン同士の Coulomb 衝突子の速度分布関数が非等方なものへと変形される可能性

によって，NB 高速イオンの速度分布関数は Maxwell 分がある．図 1（b）は非線形 Coulomb 衝突による α 粒子速度布を仮定した計算に比べて，速度方向・ピッチ角方向に分布関数の変形を示したものである．このような速度分

広がった分布を形成する（ 5 次元ドリフト運動論コード

布関数の非等方性や速度方向への正の勾配は，プラズマ

GNET による計算例は文献［14］参照）．

の不安定性を駆動することが知られており（高エネルギー

高速イオンの速度・ピッチ角が変われば粒子軌道も変粒子駆動不安定性に関してはレビュー論文［15，16］や本誌

わり，結果として高速イオンの輸送も変わるため，高速小特集［17］参照），Coulomb 衝突の非線形性がプラズマのイオン閉じ込めの正確な評価のためには Coulomb 衝突のエネルギー閉じ込め性能に影響を及ぼす可能性がある．

非線形性を考慮することが必要である．しかし，既存の3．2．3．2

中性子減衰時間に対する非線形 Coulomb 衝突の影響

磁場閉じ込め装置では，背景プラズマの密度が高速イオ

このような Coulomb 衝突の非線形性は大型ヘリカル装

ン密度に比べて十分に高い場合が多く，この非線形性が重大な影響を及ぼす場面は限られている．

置（LHD）でも実験的に観測された．ここでは，LHD のしかし，将来的に実現されるであろう核燃焼プラズマで重水素放電で行った実験を紹介する［18］．LHD の重水素

は，Coulomb 衝突の非線形性がプラズマに重大な影響を放電では，ビームイオン -- 熱イオン間の核融合反応による及ぼす可能性が考えられる．例えば，ITER の DT 実験で中性子発生率（Sn ～1014－1015 s－1）が熱イオン同士の核は，DT 核融合反応生成物である 3.5 MeV の α 粒子と，加融合反応による中性子発生率（Sn ～1011－1012 s－1）に比熱・電流駆動目的で用いられる 1 MeV の重水素 NBI 由来

べて十分に大きい．また，ビームの入射エネルギーの範の高速重水素イオンが数十 keV の高温プラズマ中に共存

囲（Einj≲180 keV）において，DD 核融合反応の反応断面することになる．Coulomb 衝突の非線形性を考慮しない

積は単調増加である．このことから，中性子発生率から場合，熱イオン同士の核融合反応で生じる α 粒子の速度分

高速重水素イオンの速度分布関数をある程度推定可能で

布関数はほぼ等方的な速度分布となる（図 1（a）粒子軌道

ある．
LHD には 3 本の接線 NBI が設置されている（図 2 ）．この 内，NB#1 と NB#3 が 反 時 計 回 り 方 向，NB#2 が 時 計
回り方向である．この実験では NB#1 と NB#2 は軽水素
NB，NB#3 が重水素 NB, プラズマの主ガス種は重水素である．図 3 に示すように，重水素 NB の入射によって，中性子発生率が増加し，定常状態に達する．重水素 NB が切れると，それに伴い中性子発生率が減衰する．この中性

図1

TASK/FP を用いたモデル計算によって得られた α 粒子の速度分布関数［18］．横軸は磁力線方向の運動量，縦軸は磁力線に垂直方向の運動量．簡単のため，磁場の不均一性による影響は考慮していない．

図2

325

LHD に設置されている NBI システムの鳥瞰図．

Journal of Plasma and Fusion Research

Vol.100, No.7

July 2024

けるこの損失時間については文献［21，22］参照）．次に，
重水素 NB と同方向の軽水素 NB を重畳した放電（丸）と逆方向の軽水素 NB を重畳した放電（四角），軽水素 NBを入射しない放電（ひし形）を比較すると，軽水素 NB を重畳した放電では，他の ２ パターンに比べて計測値 τnexp が伸びる傾向が現れている．軽水素 NB の入射による重水素高速イオンの閉じ込め時間への影響は無いと考えると，中性子減衰時間の差は軽水素高速イオンと重水素高速イ

図3オン間の Coulomb 衝突によって，重水素イオン速度分布

典型的な放電波形．（a）は NB のポート出力，（b）は中性子発生率 . この放電は重水素 NB と同方向の軽水素 NB を重畳した例．

関数が速度方向へ拡散するためであると考えられる．図 5は 3 次元ドリフト運動論コード TASK/FP による計算結果
で，重水素 NB 方向のピッチ角（θ ～0.17 radian）に対して，子発生率の減衰は主に二つの要因に分けられる．
一つは重水素 NB 高速イオンの減速による核融合反応断

高速重水素イオンの速度分布関数を表示したものである．

面積の低下，もう一つは何らかの要因による NB 高速イオ

軽水素 NB 高速イオンとの非線形 Coulomb 衝突によって，ンの損失である．ビームの減速に起因する中性子発生率

重水素 NB 高速イオンの速度分布関数が，より高エネルの減衰には解析的な時定数が与えられている［19，20］．こギー側にテイルを伸ばしていることがわかる．この結果，の時定数 τncl は以下の式で表される．

TASK/FP によるシミュレーション上でも，中性子減衰時（）

間に有意な差が現れる．

E03/2＋EC3/2    ，
E13/2＋EC3/2（3．2．19）

3mfastTe3/2τse≡
4 2π nee4me1/2 ln Λ

（3．2．20）

τse
ln
τncl＝
3

謝

辞

本章執筆には科学研究費助成事業 18K03587 の助成も受けました．

ここで，τse は高速イオンの電子衝突による減速時間，E0

は NB の入射エネルギー，E1 は DD 核融合反応断面積が E0

参考文献に対して 1/e になるエネルギー，EC は高速イオンの減速が［ 1 ］M.N. Rosenbluth et al., Phys. Rev. 107, 1 (1957).
［ 2 ］A. Mollén et al., Phys. Plasmas 22, 112508 (2015).
［ 3 ］S.I. Braginskii, Sov. Phys. JETP 6, 358 (1958).
［ 4 ］S.P. Hirshman and D.J. Sigmar, Nucl. Fusion 21,1079 (1981).
［ 5 ］C.T. Hsu et al., Phys. Fluids B 4, 4023 (1992).
［ 6 ］K. Ida, Plasma Phys. Control. Fusion 40, 1429
(1998).
［ 7 ］K. Nishioka et al., Phys. Plasmas 23, 032511 (2016).
［ 8 ］M. Nunami et al., Plasma Fusion Res. 12, 1203039(2017).

電子衝突による寄与とイオン衝突による寄与が等しくなるエネルギー（いわゆる臨界エネルギー），ln Λ は電子 イオン間の Coulomb 対数である．ここでは純重水素プラ
ズマを仮定し，EC ＝ 18.6Te とする．温度・密度の値には
規格化小半径 ρ ＝ 0.3 での値を代表値として用いる．これは，中性子発生率の空間分布のピークが大体このあたりであるからである．
重水素 NB に対し，同方向の軽水素 NB を重畳した放電，
逆方向の軽水素 NB を重畳した放電，軽水素 NB を重畳しない放電を行い，それぞれに対して密度スキャンを行った．解析解 τncl に対して中性子減衰時間の計測値 τnexp を表示したものが図 4 である．まず，軽水素 NB に依らず，全ての放電で τnexp ＝ τncl とはならない．これは重水素 NB 高
速イオンに有限の損失時間があるためである（LHD にお

図5

図4

中性子減衰時間の計測値（縦軸）と解析解（横軸）を重畳する軽水素 NB ごとに表示したもの．τncl が長い放電ほどプラズマの密度は低い．

326

TASK/FP によって計算された重水素イオンのエネルギーに対する高速重水素イオン速度分布関数．（a）は重水素
NB と同方向の軽水素 NB を重畳した放電，（b）は逆方向の放電．実線は軽水素 NB 高速イオンとの非線形衝突を考慮した場合，破線は軽水素 NB 高速イオンとの非線形衝突を考慮しない場合．時刻 t ＝ 0は重水素 NB を切った時刻を意味する．

Special Topic Article

3.2 The Coulomb Collision Operator for Problems Including Energetic Particles	［ 9 ］I.P. Shkarofsky et al., The Particle Kinetics of Plasmas (Addison-Wesley, Massachusetts, 1966), Chap.7.
［10］例えば A. Erdélyi et al., Higher Transcendental Functions Vol.1 (McGraw-Hill Book, New York, 1953),
Chap.3.
［11］S. Nishimura, Phys. Plasmas 22, 122503 (2015), 23,029901 (2016).
［12］S. Nishimura et al., Phys. Plasmas 22, 092505(2015).
［13］A. Stahl et al., Comput. Phys. Commun. 212, 269(2017).

S. Nishimura and H. Nuga

［14］K. Tahara et al., Phys. Plasmas 30, 082505 (2023).
［15］W.W. Heidbrink et al., Nucl. Fusion 34, 535 (1994).
［16］N.N. Gorelenkov et al., Nucl. Fusion 54, 125001(2014).
［17］永岡賢一：プラズマ・核融合学会誌 97, 281-285 (2021).
［18］H. Nuga et al., Nucl. Fusion 59, 016007 (2019).
［19］J. Strachan et al., Nucl. Fusion 21, 67 (1981).
［20］W. Heidbrink et al., Nucl. Fusion 28, 1897 (1988).
［21］H. Nuga et al., J. Plasma Phys. 86, 815860306 (2020).
［22］H. Nuga et al., Plasma Fusion Res. 16, 2402052(2021).

327

J. Plasma Fusion Res. Vol.100, No.7 (2024) 328-333

小特集 古くて新しい新古典輸送理論の新展開

3 ．新古典輸送計算の実際～近似モデルと新展開～

3. Neoclassical Transport Calculations in Practice:
Approximate Models and New Developments

3．3

大域的新古典計算の多粒子種プラズマへの適用

3.3 Application of Global Neoclassical Calculation to Multi-Species Plasma
藤田慶二

FUJITA Keiji
名古屋大学大学院 理学研究科（原稿受付：2024年 2 月21日）

電子やバルクのイオンとは大きく異なる物理パラメータを持つ不純物イオンの輸送に関しては，従来の新古典輸送モデルで一般的に用いられてきた近似の妥当性が破れる場合があり，大型ヘリカル装置（LHD）で観測されている未解明現象である不純物ホールの形成は，その顕著な一例であると考えられている．本節では，この現象のメカニズム解明を1つの目標として行われている，大域的効果や磁気面上のポテンシャル構造の組み込みといった，多粒子種プラズマの新古典輸送モデルおよび計算手法の拡張に関する近年の研究について紹介する．
Keywords:
impurity transport, impurity hole, stellarator, LHD, ambipolar radial electric field

3．3．1

不純物輸送

3．3．1．1

不純物の新古典輸送と不純物ホール現象の観測

核融合プラズマには，重水素イオンや三重水素イオン

不純物ホールとは，プラズマ小半径方向の座標 r に沿っといった核融合反応に直接寄与する燃料粒子（あるいはた電場（径電場 Er）が内側（小半径中心側）を向いていそれを模擬した水素イオン）の他に，装置の素材に由来る状況下で，不純物イオンがプラズマ中心部で形成する，する炭素原子や金属原子などがプラズマに混入し，イオ

窪んだ（ホローな）分布形状を指している．この特徴的ン化する．こうしたイオン種は，不純物イオンと総称さな分布は，電荷数 Z の大きな不純物ほど効果的に形成されれる．

ることがわかっている［3］．シンク項に対応するような不その呼び名の通り，不純物イオンは核融合炉の運転に純物の排出機構がない中で，中心部の分布がホローになっ

とって不要なばかりか，核融合プラズマの維持や核融合ていくということは，不純物イオンが外側方向，つまり反応の効率に対して深刻な悪影響も及ぼす．最も深刻な

径電場とは逆方向に輸送されていることを示しているが，例が，制動放射によって，プラズマのエネルギーを吐き

これは従来の素朴な理論予測に真っ向から反しており，捨ててしまうことである．その放射エネルギーは電荷数

研究者の注意をひくこととなった．

Z とともに大きくなるため，高 Z の金属イオンが混入する

新古典輸送流束は，その駆動力を用いてと，著しいエネルギー損失が起こり，プラズマが一気に（ n － T ）＋D T abΓa＝－∑ D11

冷えてしまう．一方で，そうした性質を逆手に取り，不b純物イオンのガスをプラズマの周辺部に吹き込むことで，

nb'

ZbeEr

b

b

ab Tb'
12b（3．3．1）

装置壁をプラズマの熱から守るといった活用のされ方もと表すこともできる．ここで a，b は粒子種のラベル，nb

研究されている．しかし，このように負の側面を防ぎつつ，および Tb は密度と温度で，プライム（ ' ）はその r 微分を示

不純物イオンを有益な形で利用するには，不純物イオンしている．駆動力の選択には任意性があるため，本節での振る舞いを理解し，制御できなくてはならない．とこは第 2 章での定義（2．18b）とは異なる駆動力の定義を採用

ろが，不純物イオンの振る舞いは，メインのイオンと比しており，それに対応する輸送係数を Dijab で表している．

べても理解が浅く，未だ理論的に捉えられていない側面

粒子種間の相互作用を考慮に入れるため，式（3．3．1）のも多い．そのことを示す顕著な例が，LHD で観測されて和はすべての粒子種にわたって取るが，粒子種 a 自身の電いる不純物ホール現象である［1，2］．

荷数 Za が大きければ，Za に比例する径電場の項が支配的

Nagoya University, Nagoya, AICHI 464-8601, Japanauthor’s e-mail: fujita.keiji.x6@f.mail.nagoya-u.ac.jp

328Ⓒ2024 The Japan Society of Plasma
Science and Nuclear Fusion Research

Special Topic Article

3.3 Application of Global Neoclassical Calculation to Multi-Species Plasma	

かに負であるものの，それよりも外側では正に遷移してとなり，極端な場合はaa ZaeErΓa≃D11Ta

K. Fujita

いることが示されている．一方で，背景パラメータの観測値から新古典輸送シミュレーションによって得られた（3．3．2）

径電場の分布は，プラズマ中心から端まで一貫して負のaa
はと，電荷と径電場の積に比例する形に近似できる．D11

値を取るものだった．つまり，この新古典輸送シミュレー

拡散係数で，常に正の値を取る．するとイオンの場合，ションの結果は，プラズマ中心部では確かに妥当なもの

右辺のうち Er を除く量はすべて正の値しか取らないため，の，それよりも外側の値は実験値と大きく異なっている．

粒子束の向きは径電場 Er の向きによって決定される．こ小半径の中心よりも外側の正電場の直接的影響で，径れが，不純物イオンの新古典輸送に関する素朴な予測で電場の値が負の内側領域にある不純物が外向きに駆動さある．それゆえ，径電場が内向き（Er ＜ 0）でありながら，れるということは考え難い．しかしこの電場分布のずれ

電荷の大きな不純物イオンほど外向きに輸送されることは，少なくとも計算に用いた新古典モデルにはまだ不完を示す不純物ホールの観測は驚きであったのである．

全な部分が残っており，3．1 節で言及されたような近似や

3．3．1．2

仮定の再検討が必要であることを示唆している．

乱流による輸送の可能性

新古典輸送の予測が大きく破られたことに対する 1 つの自然な反応は，新古典輸送以外の輸送を疑うことであ

3．3．2

 純物輸送解析に向けた新古典モデルの
不再考

3．3．2．1

局所近似と Φ0 近似の再検討

る． 2 章で述べられている通り，核融合プラズマ中では，
新古典輸送に加え，乱流に起因する輸送も生じる．そして，正味の流束は，新古典輸送による流束 ΓNC と乱流に起因

3．
1 節で言及された近似のうち，ここでは局所近似と，する流束 Γturb の和として
Γa＝ΓaNC＋Γaturb

磁気面上で静電ポテンシャルが一様であるという近似について検討する．

（3．3．3）

典型的な不純物ホールプラズマは高温で低密度であると表される（改めて，これら 2 つの流束の独立性についが，そのようなプラズマでは，径方向ドリフトによる大ては次章で議論されるが，当面これらの寄与を独立に計域的効果が重要になり，局所近似が破綻する傾向にある

算することが妥当であると仮定して話を進める）．そのた［7］．そして，バルクイオンの振る舞いの違いは，径電場

め，仮に新古典流束 ΓaNC が負であっても，乱流流束 Γaturbの分布の違いとして現れ，不純物の輸送にも著しい影響

がそれを凌ぐ大きさで正の値を取れば，正味の流束は正を与えうる．上述の，観測と局所近似モデルによる数値の値になる．こうしたことから，文献［4］では，炭素イオ

計算とで径電場分布が大きく異なるという事実は，大域ンの不純物ホールが観測されたプラズマパラメータに対的効果の必要性を強く示唆している．

し，線形計算の値に基づく準線形の炭素 C6 ＋ 乱流流束を静電ポテンシャルの一様性に関する仮定についても，計算したが，その値は負であった．その後，文献［5］では，磁気面上の非一様ポテンシャル，Φ1，は値としては小さ

観測値を中心に，温度や密度の勾配を広範に変更し，そくとも，高 Z の不純物イオンにとっては無視できない量にれぞれの場合について非線形のシミュレーションを行いなりうるため，再検討が必要である．すでに Φ1 の効果を

流束の値を調べた．しかしこの研究の結果も，不純物ホー含めた新古典輸送計算が行われており［8-12］，その結果，

ルが観測されるような典型的な条件下では，乱流輸送に

確かに Φ1 の有無により，不純物イオンの輸送にはっきりよる C

粒子束は内向きになることを示していた．

とした違いが出ることが確認された．ただその効果は，

径電場分布

元々内向きであった不純物の流束をさらに内向きに駆動

6＋3．3．1．3

するもので，不純物ホール現象を説明するような兆候は

文献［5］では，不安定化するのはイオン温度勾配（iontemperature gradient, ITG）モードで，装置半径 a で規見られなかった．

格化した径座標 ρ＝r/a∈［0，1］を用いると，0.4 ＜ ρ ＜ 0.7とはいえ，これらの計算はすべて局所近似モデルを用辺りであれば十分不安定になりうるという結果を得ていいて行われたものである．Φ1 の構造を決めるのは磁気面るが，この研究で対象としている放電で不純物ホールが上の電子とイオンの分布の偏りであり，主な寄与をする

最も深くなるのは ρ ＜ 0.4 の領域である．この研究がもうのはバルクイオンである．そして上述のように，不純物

1 つ示唆するのは，少なくともそれよりも外側の領域では，

ホールプラズマ中ではバルクイオンの振る舞いに大域的径電場は負ではなく正なのではないか，ということであ効果が大く影響していることが推測される．そのため，る．輸送のスケールから見て背景分布の時間変化が無視

Φ1 の効果を適切に評価するためにも大域的効果を含んだできるような準定常な状態では，新古典流束と乱流流束計算が必要であることが示唆される［13］．

を合わせた正味の流束が Γa ＝ ΓaNC ＋ Γaturb≃0 となる必要

3．3．2．2があるが，この研究の計算結果によれば径電場が負であΦ1 を含めた大域的モデル

こうした背景から，我々は大域的ドリフト運動論コードである FORTEC-3D［14］を拡張し，Φ1 の効果を含んだると，このバランスが成り立たないためである．
径電場が正であるというのは観測と矛盾するのではな

多イオン種プラズマの大域的計算を可能にした［15，16］．

かろうか？実はそうでもない．ある観測［6］では，不純物

そして，局所近似モデルを用いたいくつかの先行研究同ホールが形成される領域（ρ ＜ 0.5）での径電場の値は確様，炭素 C6 ＋の不純物ホールが観測されている水素プラズ
329

Journal of Plasma and Fusion Research

Vol.100, No.7

July 2024

マを模擬し，計算を行った．それによって得られた第一

3．3．3

結果の分析

の成果は，局所モデルによる計算と同じ背景分布を用い

3．3．3．1

物理的解釈た計算の結果，領域全体で負の値を取る径電場ではなく，

大域的コードによる計算の結果，径電場が負でも，不実験的に観測されているように，小半径に沿って負から純物ホール領域内で炭素イオンが外向きになりうること正に遷移する径電場分布が得られたことである．さらに，が確認できた．式（3．3．1）の右辺を見ると，径電場以外に径電場は負の領域でも，炭素の粒子束は外向きになるこ流束を駆動するのは，密度勾配と温度勾配である．不純とも確認された［16］．

物ホール領域では，炭素自身の密度勾配は負であるから，これらは，Φ1 の効果を含めていない場合でも得られた

径電場同様に流束を内向きに駆動する．また，他のイオ結果であるが，Φ1 を入れた場合，外向きの炭素の粒子束

ン種の密度はほぼ平坦であるため，不純物の流束にはほはより強く駆動されることがわかった．一方，背景パラ

とんど寄与しない．すると，径電場（および炭素の密度メータをいくつか人為的に変更し，炭素の粒子束が内向勾配）の駆動に逆らって，外向きに炭素を駆動しうるのは，きになった場合を調べると，Φ1 は炭素流束をより内向き

温度勾配しかない．実際に，不純物ホールプラズマは中に駆動することが確かめられた［17］．

心部での大きな温度勾配に特徴づけられるもので，このまた，この大域的新古典計算の結果は，先述した乱流計算結果は，その駆動力としての影響が，径電場の影響

輸送計算の結果［5］と比較しても，妥当性の高いものであに匹敵することを示している．そして，わずかでも外向ることが示された．つまり，水素や炭素イオンの新古典

きの流束を実現できれば，静電ポテンシャルの磁気面上流束と乱流流束が同オーダーかつ逆負号であることによでの非一様性が流束をより強く外側に駆動する．これにり，背景密度分布の変化が，計算対象としている現象よより，不純物ホールの形成は，新古典輸送理論の枠組み

りも十分ゆっくりであるという仮定に整合する．

内で説明できる見込みが高まってきた．

Φ1 が不純物イオンを駆動する方向は，主に水素の密度

しかし，上記の結果は，Z ＝ 6 の炭素イオンについて得分布の磁気面上での偏りで決まる Φ1 の位相と，炭素密度られた結果であることに注意しないといけない．核融合の磁気面上の分布の重なり合いによって決まる．そのた

プラズマ中には，鉄やタングステンといった，電荷がよめ，Φ1 が不純物流束の方向を変えず，大きさを増やすとり大きな不純物イオンが含まれる．そして，実験的にはいう性質は，この計算結果に限らない一般的なものとい

電荷の大きな不純物イオンほど，不純物ホールを効果的える［17，18］．よって，Φ1 それ自体では，不純物ホールのに形成すること，すなわちより強く外側に駆動されるこ形成を説明できそうにはない．しかし，Φ1 を含めた計算とが確認されている．ところが，径電場の影響は電荷 Z にの方が，乱流による流束との釣り合いの精度が高く，粒比例することは先述べたとおりである．対して，温度勾子束の定量的な見積もりの精度を上げるには，Φ1 を含め配の影響が電荷に比例して大きくなるという一般的な性る必要があることが計算により示唆された（図 1 ）．また，質はない．したがって，拡張したシミュレーションモデ

電荷が大きいほど Φ1 の影響を受けやすいため，高 Z の不ルによる計算で，炭素イオンでは径電場の影響と温度勾純物ほど効果的に不純物ホールを形成するという性質に配の影響が同程度になったとしても，電荷がより大きなは，この Φ1 による影響が関与している可能性もある．

不純物イオンでは結局径電場の影響が勝ってしまい，実験とは矛盾する結果が得られるように思える．
3．3．3．2

高衝突領域の不純物輸送

そこで考えられるのが，衝突頻度の影響である．新古典輸送を流体的なレベルで分解した式（2．22）を思い出すと，径電場の影響は，装置の大円方向の対称性の欠如（非軸対称性）を反映する粘性項に由来するのであった．しかし，衝突頻度が増してくると，衝突による摩擦項の影
響が支配的になって，軸対称な装置における状況に近づく．そして，いくつかの条件が満たされる理想的な場合には，径電場の影響が完全に消失する．これは，式（3．3．1）
において，径電場の係数がゼロになること，すなわちZbab∑ T D11 ＝0
b

図1

b

（3．3．4）となることに対応する．これは通常，電子の影響が無視

異なるコードによる炭素流束の計算値（新古典流束は逆符号）の比較図．左から局所的新古典コード PENTA による
新古典流束，大域的新古典コード FORTEC-3D（F3D）による新古典流束（白抜きが Φ1を含めない場合，塗りつぶしが Φ1を含めた場合の値），そしてジャイロ運動論的コー
ド GKV による，乱流流束である．なお，乱流流束の値は［5］の結果に，フラックスマッチングと呼ばれる補正を施した値である．

できる場合に起こりうる．また，イオン種はすべて同じ温度を持つと想定できるため，ある不純物イオン種 z に対して，径電場の係数が消失する条件はzI∑ ZI D11 ＝0
I

330（3．3．5）

Special Topic Article

3.3 Application of Global Neoclassical Calculation to Multi-Species Plasma	K. Fujitaとなる．添え字 I はイオン種を意味している．加えて，上は容易に克服できた．というのも，径電場さえ与えれば，に述べた理想的な状況とは，バルクイオンが 1/ν 領域にあ輸送係数が依存するのは温度や密度といった背景パラ

り，かつ装置が最適化配位と呼ばれる新古典輸送を抑制

メータそのものの値であって，それらの勾配には依存しするよう設計された装置であるか，磁気ドリフトに比べない．そのため，温度や密度などのパラメータの値をそて E×B ドリフトが十分大きいこと，そしてまた，Φ1 の振のままに，特定のパラメータの勾配を 0 に設定して計算幅が非常に小さいといった条件が満たされる状況である

すれば，その勾配に対応する駆動力を除いた場合の流束［19-23］．

が計算できる．これを繰り返すことで，任意の輸送係数

文献［22］は，問題の条件下における不純物輸送の解析を求めることができる．しかし，プラズマ領域全体を同的な表現を導き，Φ1 の存在（に対する不純物イオン密度時に計算する大域的モデルでは，このようなパラメータの応答）が，径電場の係数のキャンセルにどう影響を与の操作は許されない．それに，局所近似モデルと比べ大えるか調べている．その結果，わずか（eΦ1/TI ～5×10－3）

きな計算コストを要求する大域的モデルでは， 1 ケース毎でも Φ1 が含められることで係数のキャンセルは阻まれ，に複数回の計算が必要となる手法は効率の観点から見て負の径電場により高 Z なイオンほど内向きに強く駆動さも採用しづらい．

れてしまうというものであった．ただし，この計算ではそこで我々は，大域的モデルにも適用できる計算法を

Φ1 がなく，径電場の係数が完全にキャンセルする場合で新たに考案してコードに組み込んだ［24，25］．その方法は，も，不純物の流束が内向きになるパラメータ領域の計算

流束の 2 つの表現式（3．3．1）と式（3．3．6）を解析的につなであったことや，Φ1 を含めた場合でも Z ＝ 24 のケース

ぐ表現を採用することだ．新古典輸送は，時間変化の意

で，わずかながら外向きになるパラメータの組み合わせ

味でも，空間的非一様性の度合いの意味でも，比較的穏があったこともあり，この結果からは，Φ1 が含まれれば

やか背景パラメータの勾配に駆動される輸送を扱う．こ径電場の係数はキャンセルせず，必ず不純物流束は内向

れはミクロな観点から見ると，分布関数の局所平衡からきになる，といった単純な結論は得られない（先述の通のわずかなずれがり，Φ1 は流束の向きを変えず大きさを強める傾向にある）．
また解析的な表現を導く際に不純物の自己衝突の影響をδfa＝∑ ∑ δgab, j Xb, jb

無視していたり，Φ1 を局所コードで求めていることなど，
大域的コードで自己無撞着に計算を行った場合に結果が（3．3．7）jのように書けることを示唆している．ここで変わってくる見込みは多分に残されている．
いずれにしても，現実のプラズマでは，これらの条件

Xb,1＝

がすべて完全には満たされることはないから式（3．
3．4）の右辺が厳密に 0 にはなることはほぼないであろう．しか

nb'
ZbeEr
Tb'
－
，Xb,2＝nb
Tb
Tb

は駆動力で，δgab, j はそれらに対する応答を表している．

しそれでも，部分的なキャンセルであれば期待はできる．

展開式（3．3．7）を式（3．3．6）に入れれば，式（3．3．1）が得

通常であれば径電場の影響がより大きくなると考えられられる．このとき輸送係数はる高 Z の不純物ほど，その係数が小さくなりうるという〈d vx v ・∇rδg 〉

ab
Djk
＝－

この効果が不純物ホールプラズマ中で働いているならば，
不純物ホール現象の電荷依存性を説明する助けになるか

3

j－1
da
a

ab, k

/
によって計算される．ここで xa ＝ mav2（2T．
a）

もしれない．

従来の計算手法では，式（3．3．7）の左辺に対応する量を

3．3．4

電荷依存性の分析

3．3．4．1

輸送係数の評価

計算していたが，右辺の各成分 δgab, j Xb, j を計算するようにすることで，大域的コードである FORTEC-3D でも，輸送係数が評価できるようになった．

上で述べた，高衝突領域にある不純物イオンに対して，径電場の正味の係数が小さくなるという性質が不純

ただし，δgab, j も Er 依存するため，式（3．3．7）は見かけ

物ホール現象に関与しているのか．この可能性を検証す

通り駆動力に対する厳密に線形な応答とはなっていない

るためには，輸送係数を評価する必要がある．しかし，ことには注意しないといけない．また，大域的モデルは

FORTEC-3D は，式（3．3．1）の右辺のように，駆動力の線

磁気面と垂直な方向への粒子移動を考慮に入れているた形和の形ではなく，左辺の流束の値をめ，厳密にいえば，ドリフト運動論方程式にはそれに対〈d vv ・∇rδf 〉

Γa＝

3

da

a

応する項が含まれ，その結果プラズマの密度や温度の径方向の分布に変化が起こりうる．しかし，通常の古典輸（3．3．6）

送計算の対象は，背景分布の勾配が緩やかなプラズマにのように分布関数を積分して求める手法に基づいている．

おける，背景パラメータの変化よりも十分速いスケール

ここで，vda は式（2．6）で定義されるドリフト速度である．

での現象である．大域的効果はむしろ，粒子の捕捉，非そのため，それぞれの駆動力がどのように流束に寄与

捕捉の振る舞いを正しく評価するための役割が大きい．
そのため，背景パラメータの勾配が十分緩やかであれば，しているのかを直接評価して確かめることができない．

局所近似と同様に，式（3．3．7）のような展開が成立する．

同じ数値手法でも，局所近似モデルであれば，この点
331

Journal of Plasma and Fusion Research

3．3．4．2

衝突演算子の改良

Vol.100, No.7

July 2024

る不純物はなく，ホローな分布が維持されやすくなる．

これにより，大域的効果を含んだ輸送係数の計算が可

しかし，これまで大域的モデルで計算できた不純物イ

能になった．しかし，それだけではまだ十分ではない．

オンは（バルクイオンに近い振る舞いをする Z ＝ 2 のヘリ

高衝突領域の輸送を正しく取り扱うには，衝突項の改良

ウムを除けば）Z ＝ 6 の炭素イオンのみである．それ以上も必要である．現在 FORTEC-3D で採用されている衝突

Z が大きく高衝突領域に属する高 Z イオンを含めた計算を演算子は，高 Z 不純物が属するような高衝衝突領域の計算行いこのシナリオを検証するためには，衝突項の改良がに適していないためである［26］．この点を克服し，電荷必要である．また，次章で議論されるような乱流との相数 Z や衝突領域に対する径電場係数の分布のようなもの

互作用が不純物ホール形成に関与している可能性も排除が計算できるようになれば，不純物ホール現象の電荷依できていない．他にも，中性ビーム入射（NBI）が誘導す
るポテンシャルの非一様性［27］や，NBI によるトルクの存性に関する上述の仮説がようやく検証可能になる．
その結果，実際に径電場の係数のキャンセルにより，効果［28，29］，慣性力など，不純物イオンの輸送に無視で高 Z の不純物イオンの方が効果的な掃き出しが起こるこきない影響を与える可能性はあるが，十分検証されていと（少なくとも［3］によって実験的に示されたように，炭ない要因はいくつも残されている．

素 C6 ＋よりもネオン Ne10 ＋の方が効果的に掃き出されること）が大域的数値計算で再現できれば，不純物ホール現

参考文献

象のメカニズムの大部分を説明できたといえるようにな［ 1 ］K. Ida et al., Phys. Plasmas 16, 056111 (2009).
［ 2 ］M. Yoshinuma et al., Nucl. Fusion 49, 062002 (2009).
［ 3 ］M. Yoshinuma et al., NIFS Report NIFS-965, National Inst. for Fusion Science (2010).
［ 4 ］D.R. Mikkelsen et al., Phys. Plasmas 21, 082302(2014).
［ 5 ］M. Nunami et al., Phys. Plasmas 27, 052501 (2020).
［ 6 ］T. Ido et al., Plasma Phys. Control. Fusion 52, 124025(2010).
［ 7 ］B. Huang et al., Phys. Plasmas 24, 022503 (2017).
［ 8 ］J.M. García-Regaña et al., Plasma Phys. Control. Fusion 55, 074008 (2013).
［ 9 ］J.M. García-Regaña et al., Nucl. Fusion 57, 056004(2017).
［10］J.M. García-Regaña et al., Plasma Phys. Control. Fusion 60, 104002 (2018).
［11］A. Mollén et al., Plasma Phys. Control. Fusion 60,
084001 (2018).
［12］J.L. Velasco et al., Plasma Phys. Control. Fusion 60,
074004 (2018).
［13］K. Fujita et al., Plasma Fusion Res. 14, 3403102(2019).
［14］S. Satake et al., Nucl. Fusion, 45, 1362. (2005).
［15］S. Satake et al., Comput. Phys. Commun. 255,
107249 (2020).
［16］K. Fujita et al., Nucl. Fusion 61, 086025 (2021).
［17］K. Fujita et al., J. Plasma Phys. 86, 905860319(2020).
［18］I. Calvo et al., J. Plasma Phys. 84, 905840407 (2018).
［19］P. Helander et al., Phys. Rev. Lett. 118, 155002(2017).
［20］S.L. Newton et al., J. Plasma Phys. 83, 905830505(2017).
［21］S. Buller et al., J. Plasma Phys. 84, 905840409(2018).
［22］I. Calvo et al., Nucl. Fusion 58, 124005 (2018).
［23］M.F. Martin and M. Landreman. J. Plasma Phys. 86,
905860317 (2020).
［24］K. Fujita and S. Satake. Plasma Fusion Res. 17,1403065 (2022).
［25］K. Fujita and S. Satake. Phys. Plasmas 29, 123903(2022).

るかもしれない．
しかし，最後に簡単に触れるように，不純物イオンの輸送に潜在的な影響を与える効果を持つ因子は他にも存在しており，高衝突領域における輸送係数の評価が，別のシナリオを模索する必要性を示す可能性も，まだ十分残っている．

3．3．5

本小節のまとめ

複数イオン種計算用に拡張された大域的モデルの応用
例として，LHD で観測されている不純物ホールプラズマの解析を行った．その結果，実験と質的に一致する小半
径に沿って負から正に径電場分布と，不純物ホールが形成される領域で外向きであり，かつ乱流計算とも整合す
るオーダーの炭素流束が得られた．また，静電ポテンシャ
ルの磁気面上の非一様成分 Φ1 は乱流流束との粒子バラン
スを向上させた．これらは局所近似モデルでは得られなかった成果であり，ある条件下では，新古典輸送において大域的効果が決定的に重要になるということを，実験結果や異なる輸送過程との比較を以て明確に示すものとなった．
これらの結果のみでは不純物ホール現象のメカニズムを完全に明かしたとは言えないが，今後の研究の指針を与えるものとして，次のようなシナリオが浮かび上がってきた．径電場による駆動力は，ある程度の大きさの電荷を持つ不純物イオンには超えられるはずのない障壁であるというのが，従来の新古典モデルから得られる素朴な印象であった．しかし，衝突頻度の十分高い不純物にとってその寄与は，衝突の効果によって径電場の係数が小さくなることにより，条件によっては無視できるほど小さなものになりうる．そして，ひとたび大きな温度勾配による駆動が，正味の径電場の駆動を凌げば，磁気面上の非一様ポテンシャル Φ1 の効果で，電荷の大きな不純物イオンほど強く外向きに吐き出され，深い不純物ホールが形成される．中心部を離れれば電場は正であるから，その領域で径電場の係数がわずかでも残れば（完全にキャンセルすることは考え難い），それ以上内側に侵入してく
332

Special Topic Article

3.3 Application of Global Neoclassical Calculation to Multi-Species Plasma	［26］H. Sugama et al., Phys. Plasmas 26, 102108 (2019).
［27］H. Yamaguchi and S. Murakami. Nucl. Fusion 58,
016029 (2017).
［28］Y. Nakamura et al., Nucl. Fusion 57, 056003 (2017).

K. Fujita［29］M. Nunami et al., Plasma Fusion Res. 12, 1203039(2017).

333

J. Plasma Fusion Res. Vol.100, No.7 (2024) 334-338

小特集 古くて新しい新古典輸送理論の新展開

4．新古典輸送と乱流輸送の相互作用
4. Interaction between Neoclassical and Turbulent Transport
今寺賢志

IMADERA Kenji
京都大学大学院 エネルギー科学研究科（原稿受付：2023年12月 7 日）

これまで新古典輸送と乱流輸送はスケール分離を仮定してそれぞれ個別に取り扱うことが多かったが，Larmor 半径程度のミクロスケールから小半径程度のマクロスケールまでを第一原理的に取り扱える大域的ジャイロ運動論モデルが進展するとともに，両者の相互作用が明らかとなってきている．本章ではそのような新古典
輸送と乱流輸送の相互作用について，現在進行形で進められている研究の最前線の話題についていくつか紹介する．
Keywords:
neoclassical transport, turbulent transport, global gyrokinetic simulation, particle transport, momentum transport

4．1

緒言

これらの独立した計算によって得られた新古典輸送と磁場核融合プラズマの閉じ込め性能は，本小特集の本

乱流輸送の総和によって輸送を評価する方法は，スケー

題である新古典輸送に加え，微視的不安定性に由来し

ル分離を利用した妥当なアプローチであると考えられてた乱流輸送（異常輸送）によっても大きく左右される．

きたが，近年，両者の相互作用が重要となる物理現象が一般に両者は，ジャイロ運動に伴う平均化効果（有限いくつか報告されている．例えばその 1 つが，両極性条Larmor 半径効果）を保持しつつ，その旋回中心座標にお件によって決定される径電場形成である．スケール分離

した解析では，それぞれ独立した両極性条件 ∑sqsΓsneo ＝ 0，

ける速度分布関数の時間発展を，実空間で大域的に追跡∑sqsΓsturb ＝ 0（qs：粒子種 s の電荷，Γsneo：粒子種 s の新古するいわゆる「大域ジャイロ運動論モデル」［1］によって統一的に記述できる．しかしながら，乱流輸送現象の時

典粒子束，Γsturb：粒子種 s の乱流粒子束）を満たすが，ひ空間スケールは新古典輸送現象に比べて短いため，両現とたび大域的ジャイロ運動論モデルに立ち返ると，旋回

象を直接同時に取り扱うシミュレーションを実行するた

中心座標における両極性条件はめには膨大な計算資源が必要となる［2］．そこでこれまでは，δ＝ρi/L（ρi：イオンジャイロ半径，L：背景分布の－∑s

径方向の特性長）に関して 1 次のオーダーの粒子種 s の分

ρts2 ∂Er
neo
turb＝4π∑ q（Γ
）
s
s ＋Γs
2
λDs
∂t
s

（4．1）

布関数 fs1 を，ランダムな乱流揺動に関して統計平均（アとなり，両者は径電場 Er の形成を介して相互作用する可ンサンブル平均）した成分〈 fs1〉ens とそのまわりの揺動部

能性がある（λDs：粒子種 s のデバイ長，ρts：熱速度で与fs1 を用いて fs1＝〈fs1〉ens＋fs1 と分解して，それぞれの時
分

えられる粒子種 s のジャイロ半径）．また，形成される径間発展を追跡することで，新古典輸送と乱流輸送を個別

電場は，帯状流と同様に微視的乱流に対して作用するこに評価してきた［3］．結果として前者は，電磁場の時間変とから，再帰的な相互作用も無視できない．

化が δ に関して 2 次の微小量として評価されていることか

そこで近年，上述のように分布関数をスケール分離すら，磁場ドリフトに起因した定常的な新古典輸送成分をることなく，「大域的ジャイロ運動論モデル」に沿って第与える「ドリフト運動論モデル」に帰着される．一方後一原理的に解析するいわゆる全分布関数（full-f ）ジャイ

者は更に，磁力線を横切る方向の乱流揺動の相関長が巨ロ運動論シミュレーションを行うことで，新古典輸送と視的なプラズマのサイズに比べて十分小さいと仮定する

乱流輸送を自己無撞着に取り扱う試みが行われている．

ことで，閉じ込め磁場や密度，温度などの平衡分布の値

本章ではそのような両者の相互作用に関して，full-f ジャとそれらの勾配を一定とする局所近似を用いた「局所ジャイロ運動論シミュレーションによって現在進行形で行わ

イロ運動論モデル」に帰着されることが多く，この場合，れている研究の最前線についていくつか紹介する．

磁場ドリフトによる輸送が 0 となり，E × B ドリフトによる非定常的な乱流輸送成分のみが与えられる．
Graduate School of Energy Science, Kyoto University, Uji, KYOTO 611-0011, Japanauthor’s e-mail: imadera.kenji.7z@kyoto-u.ac.jp

334Ⓒ2024 The Japan Society of Plasma
Science and Nuclear Fusion Research

Special Topic Article

4．2

4. Interaction between Neoclassical and Turbulent Transport	

K. Imadera 純物輸送における新古典輸送と乱流輸送
不の相互作用

結果，非両極性径電場を介した乱流輸送と新古典輸送の

高価数の不純物に対しては新古典粒子輸送が重要であて新古典輸送が増幅される物理機構を示している．まずることが知られているが，その過程において乱流輸送とイオン温度勾配（ITG）乱流による突発的な粒子輸送がの相互作用が果たす役割について解析するため，フラン

イオンと電子で異なることによって，非両極性が生じる

スの CEA のグループによって開発された full-f ジャイロ

（図 1（a））．その結果，図 1（b）に示すように（m，n）＝（0，0）

運動論コード GYSELA［4］によって，バルクイオンとしての径電場 Er が形成され，ポロイダル方向に E × B ドリフ

相互作用が明らかとなった［8］．図 1 は，乱流輸送によっ

ト vE,θ＝－Er/B＝－（1 ＋ r cos θ/R0）Er をもたらす . ここで，重水素，トレーサー不純物としてヘリウム，ネオン，またはタングステンを用いた断熱電子を仮定した full-f ジャ

r と R0 はそれぞれプラズマの小半径と大半径を表す．この

イロ運動論シミュレーションが行われた
［5］
．特に，
（A）
非

E×B ドリフトがポロイダル方向に圧縮性を有しているこ

軸対称成分を各タイムステップで 0 とすることで新古典輸とから，図 1（c）に示す（m，n）＝（1，0）の密度揺動が過渡的送のみを考慮したシミュレーション，（B）異粒子間衝突をに励起され，その位相が磁場ドリフトと一致することで，0 とすることで乱流輸送のみを考慮したシミュレーション，

磁場ドリフトに起因した新古典輸送が増幅されることが

（C）
（A）
や
（B）
のようなフィルタリングを用いない包括的な

明らかとなった．また，その振幅は質量が大きい不純物

シミュレーションを行うことで，
（A）
と
（B）
の不純物輸送のほど大きく，向きは乱流粒子輸送におけるイオンと電子和が（C）と一致するかを解析した．その結果，半径位置にの粒子束の違いを相殺するように決まることが示された．

よっては 2 倍のオーダーで不一致が見られた．これは乱流したがって，タングステンのような高 Z 不純物を取り扱う

による Reynolds 応力に起因した（m，n）＝（1，0）のモードが場合，このような相乗効果が本質的な役割を果たすこと原因であると結論付けられており（m と n はそれぞれポロが予想される．

イダル，およびトロイダル方向のモード数），新古典輸送

4．3と乱流輸送の相互作用の重要性が full-f ジャイロ運動論シ
ミュレーションで初めて報告された．

 古典輸送と乱流輸送の相互作用によるバ
新ルクイオンの粒子ピンチ

また井戸村氏によって開発された full-f ジャイロ運動論

また，京都大学において開発された full-f ジャイロ運動

コード GT5D［6］によって，バルクイオンとして重水素，

論コード GKNET［9］によって，前節で述べた乱流輸送とトレーサー不純物としてヘリウム，ベリリウム，炭素，

新古典輸送の相互作用と同様の物理メカニズムで，バルまたはアルゴンを用いたハイブリッド電子［7］を仮定した

クイオン密度分布のピーキングが起きることが示された

full-f ジャイロ運動論シミュレーションが行われた．その［10］．図 2 は，軽水素とハイブリッド電子を用いた full-f

図1

乱流輸送によって新古典輸送が増幅される物理機構．（a）突発的な乱流粒子輸送が駆動された際の乱流構造，（b）その際に形成される非両極性径電場，（c）アルゴン不純物の摂動密度分布をそれぞれ示している［8］．

図2軽水素とハイブリッド電子を用いた full-f ジャイロ運動論シミュレーションにおいて，初期温度勾配を（A）（R0/LTi，R0/LTe）＝（10，10）
としてイオン / 電子加熱を行った場合と，（B）（R0/LTi，R0/LTe）＝（10，4）としてイオン加熱のみを行った場合に得られたイオン密度分布．
熱ソースとエネルギーシンクの半径分布がそれぞれの図内に，対応する線形 δf ジャイロ運動論シミュレーションで得られた n ＝ 15の固有関数が右側に示されている．

335

Journal of Plasma and Fusion Research

Vol.100, No.7

July 2024

Γiturb－Γeturb（0）を相殺するように，つまり粒子束の総和が 0 となるように（A）
では内向き，（B）
では外向きになっていることがわかる．加えて，n ＝ 0 の E×B ドリフトによる粒子輸送もまた，粒子束の総和が 0 となるように決定されていることが新たにわかった．これは（m，n）＝（1，0）の対流セルが主に担っており，磁場ドリフトによる粒子輸送
は非線形飽和後に 0 に漸近するが，n ＝ 0 の E×B ドリフ
トによる粒子輸送は飽和後も有限のレベルを維持し，継続的な密度変化に影響を与えていると考えられる．したがって，イオン加熱によってイオンの乱流粒子ピンチを直接的に駆動することに加え，電子加熱によって電子の乱流粒子ピンチを駆動することで，両極性条件を介して間接的にイオンの新古典粒子ピンチをもたらし，相乗的な密度ピーキングが達成できることが明らかとなった．
この現象もまた，まさに乱流輸送と新古典輸送の相互作用が重要な役割を果たしている一例であると言える．
図 3 （A-1）（B-1）（A）と（B）における n ≠ 0の E×B ドリフトによるイオン / 電子乱流粒子束の時間発展と，（A-2）
（B-2）
（A）
と（B）における n ≠ 0の E×B ドリフトによるイオン／電子
乱流粒子束の総和と，磁場ドリフトと n ＝ 0の E×B ドリ
フトによる新古典粒子束の時間発展．

4．4

新古典径電場と乱流輸送の相互作用

新古典輸送と乱流輸送が相互作用する別のメカニズム
として，径電場が乱流輸送に与える影響が挙げられる．
一般に径電場は，衝突の時間スケールで力学的平衡

ジャイロ運動論シミュレーションにおいて，初期温度勾Er－vφBθ＋vθBφ＝配を（A）（R0/LTi，R0/LTe）＝（10，10）としてイオン / 電子加熱を行った場合と，（B）（R0/LTi，R0/LTe）＝（10，4）としてイオ

1 ∂p

nq ∂r

（4.2）

を介して決定される．ここで，q，n，p，Er はそれぞれ

ン加熱のみを行った場合に得られたイオン密度分布を示

電荷，密度，圧力，径電場，vθ，vφ，Bθ，Bφ はそれぞれしている．（B）では密度分布が外向きの粒子束によって緩ポロイダル流速，トロイダル流速，ポロイダル磁場，ト

和している一方，（A）では密度ピーキングが起きているこ

ロイダル磁場である．軸対称トカマクのように新古典輸とがわかる．また図 1 に示されている固有関数から，両送が両極性を有する場合（2．4 章），ポロイダル流速はケースにおいて正のバルーニング角を持つ ITG モードがvθ ＝ k（v＊）∂rT で与えられることが知られており，k（v＊）は励起されており，線形成長率も同レベルであることが明新古典輸送理論によって決定されることから，本章ではらかとなった．これは，モードの種類やその不安定の強（4．2）式の関係を満たす長時間スケールの径電場を「新古さのみで粒子輸送の向きが決定していないことを示唆し

典径電場」と呼ぶ（T：温度，v＊：規格化衝突周波数）［11］．

ている．

局所ジャイロ運動論モデルでは一般にこの新古典径電場

図 3（A-1）と（B-1）は，（A）と（B）における n0 の E×Bが考慮されていない一方，ドリフト運動論モデルでは乱ドリフトによる乱流粒子束の時間発展を，イオン成分と流を取り扱っていないことから，新古典径電場と乱流の電子成分に分けて示している．（A）では，急峻なイオン /相互作用を自己無撞着に考慮するためには full-f ジャイロ

電子温度勾配によって負の非対角熱拡散項を介して強い運動論シミュレーションを行う必要がある．

イオン／電子粒子ピンチが駆動され，その後も一定のレ

そのような相互作用が重要となる例として，メゾスケー

ベルで粒子ピンチが維持されていることがわかる．一方ルの乱流輸送現象に対する新古典径電場の役割が挙げらで（B）では，イオン粒子束は同様に若干負に振れているのれる．GT5D によって，外部熱源で駆動される ITG 乱流に対して，電子粒子束は電子温度勾配が緩やかであるたに関する full-f ジャイロ運動論シミュレーション，いわゆ

め正の対角拡散項が支配的となり，結果として外向きのる熱源駆動 ITG シミュレーションが行われ，臨界温度勾

粒子輸送になっていることがわかる．

配近傍でイオン温度分布の硬直性が発現し，熱輸送の大図 3（A-2）と（B-2）は，（A）と（B）に お け る 磁 場 ド リ フ

部分が雪崩的に径方向に伝播していることが確認された

ト と n＝0 の E×B ド リ フ ト に よ る 新 古 典 粒 子 束 の 時 間［12］．その際に形成される電場は主に力学的平衡を介し発展を示している．先述した GT5D で得られた結果［8］て形成される新古典径電場であり，結果として，雪崩的で示されているように，乱流粒子輸送がイオンと電子でな熱輸送の伝播に伴い，温度勾配と径電場も局所的な力

異なることによって形成された（m，n）＝（0，0）の径電場

学的平衡を満たしながら伝播することが明らかとなった．

Er が，（m，n）＝（1，0）の密度揺動を過渡的に励起するこ

またその伝播方向は，径電場シアの正負によって決定さとで，磁場ドリフトに起因した新古典輸送が増幅されてれることも示されており，新古典径電場が大域的な乱流いることがこの結果でも確認された．またその向きは，輸送に大きな影響を与えていることがわかった．そのよ
336

Special Topic Article

4. Interaction between Neoclassical and Turbulent Transport	

K. Imadera

部輸送障壁が形成される現象もまた，両者の相互作用が重要な役割を果たす一例であると言える．

4．5

結言

本章では，full-f ジャイロ運動論シミュレーションで新古典輸送と乱流輸送を自己無撞着に取り扱うことによって明らかとなったいくつかの物理現象について紹介した．

図 4 （A）r ＝ 0.6a0に外部から運動量を順方向，逆方向に入射した場合と入射しなかった場合の準定常状態におけるイオン
温度分布．（B）外部から運動量を順方向に入射した場合の力学的平衡に現れる各項の空間分布［9］.

まず 4．2 節と 4．3 節では，乱流による Reynolds 応力や，乱流粒子輸送によって上下ポロイダル非対称な構造が駆動された結果，新古典粒子輸送が増幅されるケースについて紹介した．なお，乱流による Reynolds 応力は新古典熱

う な 特 性 は GYSELA や GKNET に よ る 熱 源 駆 動 ITG シ輸送も増大させることが，GYSELA による熱源駆動 ITG

ミュレーションにおいても確認されており，特に熱輸送シミュレーションで明らかとなっている［19］．一方 4．4 節の雪崩的伝播の空間的な特性長を決定付ける構造として，

では，新古典径電場が雪崩的な乱流熱輸送の伝播や，運“E×B staircase”と呼ばれる局所的な輸送障壁の存在が動量輸送に影響を与えるケースについて紹介し，輸送障示されている［13，14］．

壁形成にも寄与することを述べた．

またそのような新古典径電場は，乱流による運動量以上の結果は，乱流の時間スケールで新古典輸送が過フラックスの向きを変えることが GT5D による熱源駆動

渡的に増大すること，ないし乱流の時間スケールで過渡ITG シミュレーションで明らかとなっている［15］．これ的に形成された新古典径電場が乱流輸送に影響を与える

は乱流のポロイダル対称性が新古典径電場によって変化ことが原因であると考えられる．つまり，衝突の時間スすることで，運動量フラックスの乱流成分が変化するたケールより短い乱流の時間スケールで巨視的な構造が変めである．加えて，軸対称極限では十分に小さい運動量化する場合は，スケール分離が成立せず，両者の相互作

フラックスの新古典成分が，乱流揺動によって大幅に増用を無視できないことがわかる．ただしそのような場合，

大し，準定常状態では乱流成分とバランスすることが明

実験やシミュレーションにおいて輸送を新古典成分と乱らかとなった．

流成分に切り分けることは容易ではなく，軸対称性だけ

ま た GKNET に よ る 熱 源 駆 動 ITG シ ミ ュ レ ー シ ョ ン

でない，より精緻な定義が必要である点には，注意が必要である［20］．

によって，そのような新古典径電場と乱流の相互作用を介して，内部輸送障壁が形成されることが示されている

近年の計算機性能の大幅な向上に伴い，イオンスケー［9］．図 4（A）は，熱源駆動 ITG シミュレーションにおいルと電子スケールを同時に解くマルチスケールジャイロて r ＝ 0.6a0（a0 は最外殻磁気面の小半径）に外部から運運動論シミュレーションが可能となってきたが［21］，新動量をプラズマ電流の向きに対して順方向，逆方向に入古典輸送と乱流輸送を同時に解くという意味でのマルチ射した場合と入射しなかった場合の準定常状態におけるスケール full-f ジャイロ運動論シミュレーションが，今後イオン温度分布を示している．順入射の場合のみ内部輸

新たな物理を切り開くことを期待する．

送障壁が形成され，イオン乱流熱輸送係数が新古典輸送

謝レベルまで低下することが明らかとなった．この傾向は
HL-2A 装置における弱磁気シアプラズマでの実験結果と辞

本記事の執筆にあたって，有益なコメントをいただいた

定性的に一致している［16］．また図 4（B）は，順方向入射

京都大学の石澤明宏教授，武藤幹弥氏，ならびに日本原した場合の力学的平衡における各項の空間分布を示して子力研究開発機構の井戸村泰宏氏に深く感謝いたします．

いる．この結果から，順方向の運動量入射によって力学

本稿の作成は，科学研究費［基盤研究（C）：20K03903］

的平衡を介して強い径電場シアが形成され，内部輸送障の援助を得て行われました．

壁が形成されたものと考えられる．順入射が輸送障壁形成に有効である物理メカニズムは，文献［16］で提唱され参考文献ているように，非局所バルーニング理論［17］と運動量輸［ 1 ］J. Brizard and T. S. Hahm, Rev. Mod. Phys. 79, 421(2007).
［ 2 ］X. Garbet et al., Nucl. Fusion 50, 043002 (2010).
［ 3 ］プラズマ・核融合学会編：プラズマシミュレーション（京都大学学術出版会，2018）．
［ 4 ］V. Grandgirard et al., Comput. Phys. Commun. 207,
35 (2016).
［ 5 ］D. Estève et al., Nucl. Fusion 58, 036013 (2018).
［ 6 ］Y. Idomura et al., Comput. Phys. Commun. 179, 391(2008).
［ 7 ］Y. Idomura, J. Comput. Phys. 313, 511 (2016).

送理論［18］から説明することが可能である．非局所バルー
ニング理論によると，磁気シアが正の場合，正（負）の径電場シアは負（正）のバルーニング角を与える．一方，
運動量輸送の残留応力成分は，負（正）のバルーニング角の場合に正（負）の運動量輸送をもたらす．したがって，
図 3（B）に示す様に，順回転入射の場合，急峻な運動量分布が保持されるように運動量ピンチが発生し，結果として強い径電場シアが形成されたものと考えられる．このように，新古典径電場によって乱流輸送が抑制され，内337

Journal of Plasma and Fusion Research

Vol.100, No.7

July 2024［14］W. Wang et al., Nucl. Fusion 58, 056005 (2018).
［15］Y. Idomura, Phys. Plasmas 21, 022517 (2014).
［16］D.L. Yu et al., Nucl. Fusion 56, 056033 (2016).
［17］Y. Kishimoto et al., Plasma Phys. Control. Fusion 40,
A663 (1998).
［18］Y. Camenen et al., Nucl. Fusion 51, 073039 (2011).
［19］Y. Asahi et al., Plasma Phys. Control. Fusion 61,
065015 (2019).
［20］H. Sugama et al., Phys. Plasmas 29, 052509 (2022).
［21］S. Maeyama et al., Nat. Commun. 13, 3166 (2022).

［ 8 ］Y. Idomura et al., Phys. Plasmas 28, 012501 (2021).
［ 9 ］K. Imadera and Y. Kishimoto, Plasma Phys. Control.
Fusion 65, 024003 (2023).
［10］K. Imadera et al., 29th IAEA Fusion Energy Conference, London (2023).
［11］F.L. Hinton and R. D. Hazeltine, Rev. Mod. Phys. 48,
239 (1976).
［12］Y. Idomura et al., Nucl. Fusion 49, 065029 (2009).
［13］G. Dif-Pradalier et al., Phys. Rev. E 82, 025401(2010).

338

J. Plasma Fusion Res. Vol.100, No.7 (2024) 339-340

小特集 古くて新しい新古典輸送理論の新展開

5 ．おわりに5. Summary
前山伸也

MAEYAMA Shinya
核融合科学研究所（原稿受付：2024年 2 月21日）

本小特集では，新古典輸送理論に関する最近の研究のれている場合には，非 Maxwell 平衡分布を場の粒子とし進展について取り上げた．本小特集を締めくくるに当たて取り扱う Rosenbluth ポテンシャルの計算に立ち戻る必り，各章の内容を振り返ってみよう．

要があることが述べられた．また，非 Maxwell 分布効果

第 2 章では，新古典輸送に対する概要を振り返った．

を非線形 Coulomb 衝突項として数値的に扱う手法の進展

新古典輸送理論は輸送オーダリングが成立するようなについての一例として，大型ヘリカル装置（LHD）にお

ゆっくりとした時間変化を持つ平衡分布関数の摂動［式

ける中性粒子ビーム（NBI）入射の問題を取り上げた．高（2．4）］を記述する理論であり，閉じ込め磁場のトロイ

速イオンとの非線形 Coulomb 衝突による高エネルギー側

ダル性に起因するドリフト運動の非等方性と Coulomb 衝に伸びたテイル形成によって生じると予想される中性子突による摩擦のバランスとして生じる非等方圧力テンソ

減衰時間の変化が実験的にも観測された．

ルが引き起こす輸送現象［式（2．22）］である．標準的な第 3．3 章では，大域モデルに加えて磁気面上非一様静電

新古典輸送理論で用いられる仮定として，平衡分布の空ポテンシャルが重要となる LHD における不純物ホール形

間変化は乱流揺動の波長に比べて緩やかであるとするド成の問題を取り上げた．局所モデルでは実験と一致しな

リフトオーダリングの下では，新古典輸送と乱流輸送はかった径電場分布が大域モデルでは記述できること，磁分離して記述できる［式（2．19）］．輸送だけでなく，新気面上非一様な静電ポテンシャルが高価数不純物の新古

古典トロイダル粘性によるトロイダル回転の影響，ブー典輸送を増大させ，乱流輸送との定量比較においては考

トストラップ電流による自発電流なども記述する．ま慮する必要があることを明示した．また，さらに高価数た，線形化 Coulomb 衝突項の持つ保存則や自己随伴性のに対する依存性を調べるには高衝突周波数領域に対する性質から導かれる新古典輸送の本質的両極性（intrinsic

衝突演算子の改良が必要となることも課題として挙げら

ambipolarity）と Onsager 対称性についても触れた．

れた．

第 3．1 章では，ドリフト運動論方程式の数値解法に用い

第 4 章では，従来仮定される短波長・高周波の乱流輸られる種々の近似について説明した．磁気面上の静電ポ

送と長波長・低周波の新古典輸送というスケール分離がテンシャル変位，磁場ドリフトの有限軌道幅による非局

成立しない例について，大域的ジャイロ運動論シミュレー

所効果，最低次 Maxwell 平衡からのずれによる衝突項のションに基づく報告がなされた．乱流揺動の持つ短い時

非線形効果や保存則の数値的保証，速い平均流の有無と間スケールで長波長変動が生じることで，空間スケール

いった物理モデルや近似の違いが存在し，現在も研究対の分離が成り立たなくなり，両者の相互作用が無視でき象となっている．局所・非局所計算モデルの比較の例となくなる．4．2 節では乱流から新古典への影響として，乱して，両極性径電場の解析結果を示し，磁気ドリフト軌

流揺動による密度のポロイダル方向非一様性によって新道の磁気面からのずれを抑制するような新古典最適化配古典輸送が過渡的に増大する現象が述べられた．4．3 節で位に近い場合は局所近似モデルで定量的に十分な精度のは，両極性拡散の乱流・新古典分離が成り立たない例が計算が可能であること，ただし径電場の正負が入れ替わ示された．また，4．4 節では新古典から乱流への影響としる近傍では局所近似モデルでは解が一意に定まらない一て，新古典的に形成される径電場が乱流の雪崩的伝搬や方，大域モデルでは自然に一方の解に収束するなどの違運動量輸送に与える影響が述べられた．
非専門家の一部には，「新古典輸送理論はすでに 80 年代

いがあること［図 3．1．1］が示された．
第 3．2 章では，Maxwell 平衡の仮定が成り立たない場合には体系化され，トカマクにおける理論との良好な一致がの例として加熱されたプラズマを取り上げた．ビームイ確認されるなど，標準的な理論が完成されている」とのオンなどがあり平衡速度分布が Maxwell 分布からかけ離誤解があるかもしれない．しかしながら，本小特集で取National Institute for Fusion Science, Toki, GIFU 509-5292, Japanauthor’s e-mail: maeyama.shinya@nifs.ac.jp

339Ⓒ2024 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research

Vol.100, No.7

July 2024

り上げたように，長年デファクト・スタンダードとされプラズマ”においては非 Maxwell 平衡による Coulomb 衝てきた標準理論にはそこで用いた仮定による適用限界が突の非線形効果を考慮すべき場合があることが第 3．2 章であり，非 Maxwell 平衡・有限軌道幅による非局所・磁気例示された．磁気面上非一様静電ポテンシャルは“高価

面上不均一静電ポテンシャルの寄与・乱流 - 新古典輸送相

数不純物”に対しては輸送を増大させる傾向があり，第 3．

互作用といった，従来モデルでは説明できない現象も近3 章の例では最大 2 - 3 割の輸送増大が示された．第 4 章で年報告されている．この点で，新古典輸送理論を発展さ議論された乱流が新古典的輸送に影響を与える状況の例

せていくことの必要性は論を待たないだろう．一方，現として，“高価数不純物”や“輸送障壁”がある場合など

時点での取り組みとして，非専門家が従来モデルに基づが該当する．いずれにしても，従来モデルで用いられて

く新古典輸送コードを解析ツールとして利用する場合をいる仮定が成立しているか注意する必要がある．逆に言想定すると，「結局，これまでの新古典輸送解析を置き換

えば，標準的仮定が破れているような実験条件を精査しえる必要があるのか？それはどのような場合か？」といていくことは，新古典輸送理論の研究をさらに深めていうことは関心事であろう．各章の振り返りでも述べた通くことに貢献するものであり，読者自身の新たな気づきり，非局所性は磁場ドリフト有限軌道幅の磁気面からのや研究のタネにもつながると期待される．

ずれが大きい“ 3 次元磁場配位”において影響が出やすく，最後に，本小特集をまとめるにあたり，お忙しい中執

第 3．1 章では“両極性径電場の分岐・遷移現象”の評価に筆いただいた著者の先生方に感謝し，終わりの言葉とさ違いが出る例が示された．また，“ビーム成分を持つ加熱せていただきます．

小特集執筆者紹介
さ

たけ

しん

すけ

いま

佐 竹 真 介自然科学研究機構 核融合科学研究所 構造形
成・持続性ユニット・准教授．長らく新古典
輸送シミュレーション研究に携わってきまし
たが，この数年は機械学習を用いた最適化配位探索の研究も行っています．子どもの誕生を機に筋トレに励むようになり，筋肉もだいぶつきましたがなぜかそれ以上に脂肪もついてしまい，自身の巨大化に頭を悩ます日々です．
にし むら

まえ

伸

が

ひで

お

た

けい

や

自然科学研究機構 核融合科学研究所 メタ階
層ダイナミクスユニット 准教授．プラズマ乱流の研究に関して，シミュレーション開発から物理解析，データ分析・モデリング手法に興味を持っています．子どもの小学校 PTA つながりでバレー
ボールを始めました．出張などで参加できないこともありま
すが，久しぶりに体を動かすのを楽しんでいます．

奴 賀 秀 男ふじやま しん

前 山 伸 也

核融合科学研究所 構造形成・持続性ユニット
助教 主な研究分野は，複数イオン種状態，高エネルギー粒子存在状態への新古典理論の拡張と応用．所内の様々な「研究系」を異動しながらそれをやっておりましたが，これもこのたびの体制変
更の副産物でしょうか．今回分担させていただくこととなりましたこの小特集が皆様のお役に立てば幸いです．

ぬじ

京都大学大学院 エネルギー科学研究科 プラ
ズマ・核融合基礎学分野 准教授（エネルギー博士）．専門はグローバルジャイロ運動論シ
ミュレーション．本小特集のタイトルである「古くて新しい新古典輸送」．私が考えた訳ではないですが，
良いネーミングだと思っています．私も「歳を重ねたがまだまだ若い新中堅教員」として頑張っていきたいです．

しん

西 村

でら けん

今 寺 賢 志じ

藤 田 慶 二名古屋大学 大学院理学研究科 研究員．昨年より核融合プラズマの新古典輸送から，粒子加速とオーロラの構造形成に主な研究対象が移り，プラズマ現象の奥深さと幅広さを改めて実感しています．火の玉や UFO として報告される発光現象の一部などいわゆる超常現象の中にも，プラズマ物理で説明可能と議論されているものがあるようで．ゆくゆくはそうした現象の科学的解明を目指した研究もしていけたらと思っています．

340