J. Plasma Fusion Res. Vol.90, No.6 (2014)３
３８‐３４５講座 高密度相対論プラズマの粒子シミュレーション技法

３．広密度プラズマシミュレーションのための衝突モデル
3. Collision Models for Large Density Scale Plasmas
千 徳 靖 彦， 田 口 俊 弘１）
SENTOKU Yasuhiko and TAGUCHI Toshihiro１）
ネバダ州立大学リノ校物理学科，１）摂南大学理工学部電気電子工学科（原稿受付：２
０
１
４年４月２
３日）

粒子コードにおけるクーロン衝突過程の計算手法を紹介する．レーザープラズマではレーザープラズマ相互作用領域から固体密度領域までの非常に幅広い密度領域を取り扱う．このため，密度の違いを粒子の重みにより表すことが多い．このような異なる重みをもつ粒子間の相対論的クーロン衝突過程をTakizuka-Abeモデル（２体衝突モデル）を発展させて取り扱う方法があり，一方で，２体衝突モデルとは別に Langevin 方程式に基づいて解く，統計的手法も開発されている．本節ではそれぞれの手法の紹介とともに，利点ならびに問題点などについても触れる．
Keywords:
particle-in-cell, kinetic plasma, collisional plasmas

３．
１ はじめにかし統計的な運動量とエネルギーの保存の精度を維持する

本章ではプラズマ中の衝突過程について特に広密度プラ

ためにはセルあたり１００個以上の PIC 粒子が必要となり，ズマを取り扱うための技法を紹介したい．ここで取り扱う

それ以下だと運動量／エネルギー保存の精度は悪くなり実広密度プラズマとは，例えば超高強度レーザーにより生成

用的でない．ちなみに異なる重みの粒子間の２体衝突で完された固体プラズマのようなものである．レーザーは固体

全に運動量とエネルギーを保存し，正しく緩和過程を記述

密度より低いレーザー光の臨界密度以下の無衝突プラズマするスキームは存在しない．Sentoku と Kemp は Ref.［３］の内で主に吸収され，吸収されたエネルギーは高速電子によモデルを強相対論的な領域まで拡張し，また異なる重みのり衝突過程が重要となる高密度領域に輸送され固体は加熱

粒子間衝突でエネルギーを完全に保存し，運動量を統計的される．吸収領域と固体内の密度比はターゲットの物質にに保存するスキームを開発した（２００８，SK08 モデル）［５］．

よって１００から 1,000 倍超となる．前章で議論したように粒SK08 モデルではエネルギー保存が満たされることにより

子シミュレーション（PIC）でこのような広密度プラズマ

メッシュ当りの粒子数が少なくてもシステムのエネルギーを効率的に取り扱うには重みの異なる粒子を使うことにな

保存性が高いことが最大のメリットである．次節では SKる．本章では２体衝突モデルと統計論的な衝突モデルを紹

０８モデルを紹介する．

介し広密度プラズマに適用する方法を議論する．

２体衝突モデルとは別に Langevin 方程式に基づく統計

クーロン衝突のシンプルな２体衝突モデルは Shanny 等的手法も開発されている．衝突モデルとは，結局のところ，［１］が提案し，Takizuka と Abe が複数の粒子種をもつプラ

高密度プラズマ中での速度緩和を模擬するためのものであズマに拡張した（１９７７，TA77 モデル）
［２］．TA７７モデルはり，粒子が実際にクーロン力で軌道を修正する様子を計算

非相対論的なエネルギーの衝突において，その前後で運動するわけではない．２体衝突モデルとは粒子のデータから

量とエネルギーが完全に保存する．TA７
７モデルは Sen-

直接速度緩和に相当する速度揺動を計算する手法である．

toku 等により弱い相対論的なエネルギー領域に拡張されこれは保存性が良いという利点がある反面，粒子の識別がている［３］．また Nanbu と Yonemura が２体衝突モデルを不要であるという粒子コードの特性が失われ，どの粒子が異なる重みをもつ粒子間の衝突に拡張した（１９９８，NY98

同じセル内に存在するかを探索するためのランダムペアリ

モデル）
［４］
．NY98 モデルでは重みの異なる粒子間の衝突

ングのプロセスに時間を費やす．これに対し，Langevinにおいて，後で説明するようにリジェクション法を用いて方程式を用いた場合にはあくまでも粒子の平均量である平いるので，個別の衝突では運動量とエネルギーは保存しな

均速度と平均エネルギーのみで緩和計算が可能なので，粒い．しかし何度も衝突を繰り返すうちに統計的に運度量と子コード本来の構造を失わずに計算することが可能であエネルギーが保存されるように工夫されている．NY98 はる．これについての詳細な議論は３．
３節で述べる．

等しい重みの粒子間の衝突において TA77 に一致する．しDepartment of Physics, University of Nevada, Reno, USA
1)Setsunan University, Neyagawa, OSAKA 572-8508, Japanauthors’ e-mail: sentoku@unr.edu, taguchi@ele.setsunan.ac.jp

３３８!2014 The Japan Society of Plasma
Science and Nuclear Fusion ResearchLecture Note

3. Collision Models for Large Density Scale Plasmas

３．
２ 広密度プラズマのための２体衝突モデル

Y. Sentoku and T. Taguchiてみよう．

２体衝突モデルはセル内で起る可能性がある全ての衝突

図１は重み２の粒子 #(と重み６の粒子 $)が衝突した場

過程を衝突ペアのランダム選択によって表現している．そ合である．イメージとして基準となる重み１を内包粒子１のため同種粒子間（電子‐電子など）や異粒子種間（電子‐としている．#(内の内包粒子数は相手側の $)粒子の内包粒イオン）の個々の衝突が等しい確率で選択されるように考子より少ないので，## $&"&$"となり，必ず衝突により散慮して，ランダムソーティングによるペア選択を行う．ペ

乱される．一方，$)粒子の内包粒子は確率的に６つのうちアリング手法の詳細は参考文献［２］にあるのでそちらを参２つ照のこと．また，非相対論の TA７７を相対論に拡張するに１より低い．そこで NY98 モデルでは $)粒子が衝突によりは，速度による計算を運動量ベースに変更し，衝突する２散乱される確率を 2/6 とし，３回に１回しか散乱されず，粒子の重心系への変換をガリレオ変換からローレンツ変換

他２回は #(粒子が散乱され運動量やエネルギーが変化しに変えればよい．詳細は参考文献［５］にあるのでこれも割ても $)粒子はそのまま変化しない．これが NY98 がリジェ

愛する．

クション法と呼ばれる理由である．リジェクション法はモ

しか衝突しない．#$ $#"&となり散乱される確率は

衝突周波数の扱いについて述べておく．次式は二つの粒

ンテカルロ法ではよく使われる手法であるが，この方法で

子種 #
（電荷 &
（電荷 &
#）と $
$）間の衝突における衝突周波は運動量とエネルギーは個々の衝突で保存しない．しかし

数である．

ながら多くの粒子を同じように衝突させた後にエネルギー

#
'%
&
%
#&
*!
$&,
&
．
#$ $
#
-2+/.2+/と運動量は統計的に保存する．ただし，あとに検証で示さ（１）

れるように，セル内の粒子数が１００個以上と多くないと統計上保存性が悪い．

ここで ,*は二つの電荷の密度のうち低い方の密度を示し，

さて衝突の前後でエネルギーを保存するにはどうしたら!はクーロン対数，また -2+/と .2+/は２粒子間の相対的な運いいだろうか．Sentoku と Kemp により提案された手法で動量と速度である．これはプラズマ中の衝突周波数であるは両粒子は衝突において必ず散乱される（図１（b）を参Spitzer の理論を使っている［６］．広密度プラズマでは高温

照）．ただし重い方の $)粒子は #$ $#"&の確率で部分散乱

低密度プラズマから，低温高密度プラズマを連続的に取りされるものとする．そして $)粒子は内包粒子間でエネル扱うが，低温高密度では縮退している場合 Spitzer モデル

ギーと運動量を分ける．そのため衝突後（after）の $)粒子が使えなくなる．Spitzer の衝突周波数は温度の 2/3 乗に反のエネルギーは確率 #$ を使い衝突の前（before）の値と散比例するので，エネルギーが低くなると衝突周波数は発散

乱による変化した分（scattered）の平均値をとりしてしまう．物理的にはある温度以下では縮退が起り衝突

'
,
4
+2
(+,
12
+
3
)
'
4
4
+2
+*
#$，
%
""#$&
$%
#
%
"%
#
$
$
$

周波数は発散しない．縮退状態への閾温度$!は以下の式で計算される．

（５）となり，同様に運動量は，##"$$!$ ! %
'#,'&
#
．
+'
#

'
,
4
+2
(+,
12
+
3
)
'
4
4
+2
+*
!$
#$，
"!#$&
$!$
#
%
"!$
#

（２）

（６）

となる．衝突の前後では２粒子の全エネルギーと全運動量

ここで %は光速，,'%+'&は電荷 &#と &$の密度のうち高いは保存している．しかし式（６）の衝突後の運動量を使っ方の密度（質量）を表す．この温度以下では縮退した状態

'
,
4
+2
より小
て計算したエネルギーは式（５）のエネルギー %$とし衝突周波数を以下の一定値とし発散を抑える．

さくなってしまう．つまり単純に確率#$にて運動量を分配

#
&
+'%
%
#&
$&!
&
．
#$ $
$
!
'
$

するとエネルギーを失ってしまう．ここで物理的なイメー（３）ジを考えてみると，一つのマクロ粒子で部分散乱が起ると，そのエネルギーはマクロ粒子内部の内包粒子が熱化しここから広密度プラズマを取り扱う上で大切な異なる重みをもつ粒子間の衝突について説明する．まず Nanbu 等による NY98 モデルについて説明し，その問題点を改善するための技法を紹介したい．
まず #粒子種の (番目の粒子の重みを /#-とし，同様に/$.を $粒子種の )番目の粒子の重みとする．するとグリッド内の粒子種の密度 ,# と ,$ はそれぞれ，"#,# $! /#(，("$,$ $! /$(．
)

（４）

いま #種の (番目の粒子が $種の )番目の粒子と衝突したとすると，粒子 #(の衝突確率は ## $/$)"/#(!
/$)&であ0'
5%
り，一方の粒子 $)の衝突確率は #$ $/#("
/#(!
/$)&と
0'
5%

図１

なる．わかりにくいと思うので図の例でこの確率を確認し
３３９

異なる重みをもつ粒子間の衝突．（a）
NY98：リジェクション法．（b）
SK08：マージング法．

Journal of Plasma and Fusion Research Vol.90, No.6 June 2014て内部エネルギーが増えているはずである．しかし PIC

は，&種粒子と '種粒子の数密度はそれぞれ同じである．

粒子に内部エネルギーは無いので，以下の式で重みの大き

同じ重みの（a）の場合，マクロ粒子のペア数（!',&'/.）は

い粒子 'に熱化分の運動量を便宜上足し合わせエネルギー

TA77 のルールに従い１０となり，衝突数 &&'も１０となる．異の不足分を補うこととする．

なる重みの（b）の場合マクロ粒子の衝突ペア数は４対であ

)
*
-&
+
&
)
0
(/
!'
#!'
"%!'$．

り &&' は４となる．つまり実際に起こらなければいけない（７）

衝突数１０に対してその数は 2/5 となっている．そこで時間%''$ の値はエネルギーが保存するように以下の式で決定&)
0(/
し，方向は運動量 !'に垂直なランダムな方向とする．
#
#
．
%'
'$'
#%'" )
!)
'
"(
&
)
0
(/
'

&
)
0
(/
'

（８）

ステップに%&&##&&'&$#&倍を掛けて，全衝突数を１０に補#%正する．これにより物理的な衝突数は一様な重みの場合（a）と異なる重みの場合（b）は同等になる．もし &と 'が同一種である場合，マクロ粒子のペア数 ! は !#!&##（!&

ここで %' は粒子種 'の質量であり，また，
#
&
)
0
(/
)
#
%'"#，
#""(
'
(
&
)
0
(/
'

#
&
)
0
(/ #
)
''
%'"#．
#"""'
'#
'
&
)
0
(/
'

（奇数のとき）であるのが偶数のとき）か !& #%
!&""&
#
#
で，全衝突数 &&& は

（９）

!

（１０）

&&& ##!
$%''$ をランダムな方向に設定することで，足し合わせた
時に %''$ はキャンセルし統計的に運動量保存が満たされ)#$!")#$)#$!""
)#$&，
,&
1%

（１３）

となり，時間ステップを %(#%
&&#
%(
&&&&
!とすればよい．
ここから TA77，NY98，SK08 をビーム緩和の シミュる．
NY９８には異なる重み間の衝突を確率的に正しく扱うたレーションに適用し，エネルギーと運動量の保存性についめに，マクロ粒子内の実粒子数を使った時間調整という大!
!!""の温度をもつ背景電子を考て調べてみよう．熱速度!

切なしかし判りづらい概念があるのでここで解説しておえる．密度は 1025 cm−3 とし，(#!で背景電子の内１０％にく．セル内のマクロ粒子が !& $!' であるとき，２体衝突!
%#"のドリフト運動量を与える．計算では e-e 衝突だけ!
%するペア数は !&である．また実粒子の全衝突数 &&'は以下を考慮し電磁波は解かず，空間０次元運動量３次元の純粋の式で与えられる．

なモンテカルロシミュレーションを行った．外部からエネ)&$)'$
)&$"
)'$&．
1%
$ ,&

ルギーの注入はないのでトータルエネルギーは保存しなけ!&&&' #!

ればならない．以下三つのケースを行った．

（１１）

一 様 重 み の Case A（TA77）を こ の 計 算 の 基 準 と す

全衝突数を粒子の重みが一様である場合と同等にするため

る．図３に背景電子のエネルギーの時間発展をプロットしにシミュレーションの時間ステップ %(!を以下のように補た．粒子数がセルあたり２５０個のときは NY98 も TA77 に近正する．

い結果となっている．しかし１０個になるとエラーが大きく．
%(# && %(&&' !

なっていることが確認できる．一方 SK08 は１
０個でもエ
（１
２）

&&'）とする．これの意
もし !' $!& であれば係数は（&'#味するところを図２のケースを例として説明する．図中二つのケース（a）同じ重みの衝突と（b）異なる重みの衝突で

図３

図２

サンプルケース：（a）
一様な重みをもつ場合，（b）
異なる重
み粒子間の衝突の場合．

３４０

ビーム緩和計算における背景電子のエネルギーの時間発
展．
Case A（TA77）：一様重みの粒子，５００粒子／セル（背景電子４５０，ビーム５０）．
Case B（NY98）：異な る 重 み の 粒 子，２５０粒 子／セ ル（背景電子１
２５，ビーム１２５）と１０粒子／セル（背景電子５，
ビーム５）．
Case C（SK08）：異なる重みの粒子，粒子数は Case Bと同じ．

Lecture Note

3. Collision Models for Large Density Scale Plasmas

Y. Sentoku and T. Taguchi

ラーは小さく抑えられ最終的な背景電子の温度も TA77

ホーンの広告があった．ノイズをどう抑えるか仕組みの簡と一致している．

単な説明があったが，その中にランダムノイズは同じ周波続いてシステムのエネルギーと運動量の時間発展をみて数のノイズを足し合わせることで弱められるとあった．ひ

みよう（図４）．TA77 はどちらも完全に保存する．図はセ

らめいて帰宅後，電流計算を衝突の前後で行いそれを足しルあたりの粒子１０個の時の NY98 と SK08 を比較している．

合わせて平均を行ってみたら，驚くほどエネルギー保存がSK08 は完全にエネルギーが保存し NY98 はエネルギーの改善された．試されたい．

エラーが振動し最終的に４０％近いエラーとなっている．ま３．
３ Langevin 方程式を用いた統計モデル

た SK08 は運動量でもよい保存性が確認できる．NY98 の場

ここからは２体衝突モデルとは全く別の手法である統計

合このテスト計算の場合，セルあたり５００個の粒子を使うモデルを紹介しよう．プラズマ粒子コードでは，粒子デーことで良い結果を得られた．
本節をまとめると，TA77 と NY98 を発展させた SK08

タと格子データをやり取りすることでプラズマの相互作用のエネルギー保存性と運動量の統計的な保存性を確認できを計算する．この方法の利点は，計算量の低減である．もた．前章の高次補間と組み合わせることで数値加熱や数値し粒子間力を２体相互作用で解くと，粒子数の２乗に比例

エラーを抑制し計算精度を上げることができる．これでする計算量が必要であるが，空間格子を導入して電場や磁ハッピーと言いたいところであるが，実はまだ一つ問題が場を格子点で解けば，粒子数に比例した計算量に格子点数残っている．シミュレーションは低密度領域と高密度領域に関係した計算量の加算で済む．格子点数は粒子数よりもで同じ時間ステップを使い進めるため，高密度領域では少ないので，格子点に関連した計算量が格子点数に関する$"#!!が大きくなり強く粒子が散乱されるためランダムノ

多少大きいべき乗依存性があっても問題にならない．

イズが発生し，電流にそのノイズがのってしまう．その結

ただ，格子の導入とは空間を粗視化することであり，場果電磁場と電流の一貫性が壊れ数値加熱によりエネルギーの量の短波長成分が除去されることを意味する．短波長成が増大する．これは強い散乱で発生したノイズで前章で紹分の除去は，粒子の近接的な相互作用を抑制するため，粒介した粒子数依存のノイズと異なりグリッド計算で高次補子コードが比較的大きなタイムステップでも安定に動作す

間をしても落とせない．UCLA の OSIRIS では高次補間し，る方向に働く．しかし，衝突の効果が重要な高密度プラズ

さらにグリッド間で電流の平滑化を取り入れていると聞い

マ領域では，この短波長成分の除去による近接相互作用のたが本質的な解決法ではない．このランダムノイズのため抑制は物理的に正しくない結果を生む．かといって，空間高密度のシミュレーションを行うとエネルギー保存性が極格子の間隔を狭くして短波長成分を含めるようにすると，めて 悪 く な る．時 間 ス テ ッ プ を 小 さ く し，高 密 度 域 で

タイムステップを小さくしなければならず，計算量の大幅$!!とすればもちろん抑制できるが，それでは計算量"#!!

な増加が避けられない．すなわち計算機シミュレーションが増えてしまう．これが当時の私を最後まで悩ましたが，としての実用性を損なうことになる．

解決方法はふとしたことから浮かんだ．どこかの会議の帰

しかし，プラズマ粒子コードにおいて，個々の衝突におりの飛行機で，座席のシートポケットにある雑誌をぺらぺ

ける軌道の偏向は重要ではない．そもそも式（１）のクーロらやっていたら，ノイズキャンセリング機能をもつヘッド

ン衝突係数は小角度散乱の集積結果であり，２個の実荷電粒子の衝突に要する時間は，プラズマ振動数やクーラン条件で決まる粒子コードのタイムステップに対して遙かに短く，マルコフ過程と考えて衝突計算をすることができる．
コードとしては電場磁場中での粒子運動に衝突の効果による非可逆的な速度緩和や熱緩和の効果が入ればいいのであり，これが付加的に衝突モデルを導入すれば十分であることの理由である．
前節で説明した２体衝突モデルでは，粒子データを元に衝突計算をしているが，あくまでも緩和量を模擬するのが目的なので，軌道を正確に解く必要はなく，運動量の修正だけで事足りる．このため，２体衝突で計算すべき必然性
はなく，緩和を模擬できるような他のアルゴリズムを選択することは可能である．
さて，２体衝突モデルは粒子データから緩和量を評価する方法である．そもそも粒子コードは粒子データと格子データのやり取りにより粒子間相互作用を計算する手法であるため，個々の粒子の識別は不要である．たとえば，電子密度を計算する場合，粒子の位置データを元に格子配分

図４を行うが，このとき，粒子がどの格子に属するかをあらか

全システムのエネルギーと運動量の時間発展．

じめ識別しておく必要はなく，粒子の番号順に計算するこ
３４１

Journal of Plasma and Fusion Research Vol.90, No.6 June 2014とが可能である．しかし，２体衝突モデルで計算する場合

式（１５）は，見てわかるように分布関数 &"に関しての非線は，どの粒子とどの粒子が近いかという探索を行わねばな

形方程式であり，これを直接解くのは数値解析でもかなりらず，粒子のソーティングなどの余分な処理が必要とな難しい作業である．まして，粒子コードにおける粒子デー

る．これは，粒子コードの本来の利点を損ねているといえタはセルごとの粒子数が少なく，意味のある分布関数の形

る．また，前節で詳細に説明があったように，重みの違う状を求めるのさえ困難である．

粒子が存在する場合には，これを考慮して２体衝突計算を

この係数を計算する手法はいくつか提案されているが行う必要があるが，これも粒子に個性をもたせたために生［８，９］，ここではメリーランド大学の Antonsen が提案したじた問題点であると考えられる．

簡易モデル（簡易 Langevin モデル）を紹介しよう［１０］．衝これに対し，格子データから緩和に必要な運動量の修正

突モデルは，要するに熱緩和，すなわち非平衡状態から平計算をすることを目的としたのが，Langevin方程式を用い衡状態への状態変化を解くのが目的である．今，現在の分たモデルである．Langevin 方程式とは，粒子の運動方程式

布が平衡状態からそれほど大きく外れてはいないと考えにランダムな揺動力を付加したものである．この方程式る．上記のように，この問題はあくまでも緩和係数が実際は，ブラウン運動を解析するために考えられたもので，次と大きく異ならなければよい．そこで，終状態までの熱緩式の形式である．なお，本節では簡単のため非相対論的な

和時間をある程度の誤差内で評価すればいいと考えれば，記述で話をする．

Rosenbluth ポテンシャルの計算に必要な #種粒子の分布関

&
．
&!%(
( ' %"!(%'
+
"$'
(
'+

数を，現時点ではなく終状態の平衡分布関数で評価すれば（１４）いいだろうというのが Antonsen のアイディアである．この考えに基づき，#種粒子の平均速度 &# と温度 $# で指定式（１４）で，"は電磁力などの外力であり，これに統計的平されるマクスウェル分布で Rosenbluth ポテンシャル " と均値が０の揺動力 $'&!%(を加えた+(と速度緩和項!(%'! を評価すると，速度空間フラックスは次式のようになのが Langevin 方程式である．ここで，%は速度緩和係数，

る．

%はセル内の平均速度である．揺動力は乱数を使って計算

# #
%
)#
"%
# #
#"!# %!
$#
# #
&!&#)
&$
( )
%
! "
$
( &
" '
"
&
-(
!&#(
" & # &)&
% "'
&('
．（１９）
(#
(#

するが，これが粒子散乱項になる．しかし，ランダムな力を与えるので，平均的には粒子を加速させ，時間とともに
平均エネルギーが増大する．すなわち加熱作用がある．この加熱を抑制するのが速度緩和項であり，こちらは冷却作用をもつ．この揺動力による加熱と速度緩和項による冷却

ここで，('-(は以下で与えられ，誤差関数を使って計算すが釣り合ったところで平衡分布が決まる．逆に言えば，最ることができる．

終的な緩和温度により，$'+(と %の関係が決まるように

-

%
('
-(
% ( '#%!' ''．
,- !

なっている．これはアインシュタインの関係に他ならない．
次に，緩和係数の計算について述べる．まず，運動論に（２０）

ただし，-は次式で与えられる．

よれば，衝突過程は次式で表される［７］．
"
)
&
#"!#．
%!' &!$
)
+
#

#

-%
)

（１５）

#
&!&#(
(#'
'&$# ．
#

（２
１）

ここで，&
" は速度分布関数である．"や #は粒子の種類

以下，簡単のため，電子１種類の緩和について考える．

を表す．この式は，"種の粒子の分布関数の時間変化が衝

この場合，"種と #種が同じ粒子なので，分布関数の時間

突相手の #種の粒子（標的粒子）の速度分布で決まる緩和変化は次式のように簡単化される．

量を用いて計算されることを示している．

)
& )
&
'$ )
&&
,(
% $%
" &
%*'
#
$
&
&
& ．
(( )
)
+ )

式（１５）における #"!# は，速度空間でのフラックスを与
え，
"!
#

#

（２２）

ここで，%
,(は次式で与えられる．
*'

# #
%
# "
"%
"
#
"
&(
&(
&(
&(
&'
&
! &)$
+
%
&)"'
*
'
&)&)!'
&，
# #%
#
(
&$
%
"
!

,(
%
%
*'

（１６）
である．ここで，" や ! は Rosenbluth ポテンシャルと呼

)(%%#
-(
'
．
# # $(
&$
(,
%
! (

（２３）

この分布関数の時間変化に関する方程式 は，FokkerPlanck 方程式，

ばれ，次式で与えられる．
#
&-( &
'
(
&(
"'
'&
%!"
" ""( (#&!&) ，
)

（１７）

#
&-(
&-!&)
&(
!'
'&-，
'
)
%(&

（１８）

)
&
) "& ) ) !
%! $
' (
"
(
：' &
，
&
&)
&
)
)
)
+

（２４）

と比較することで，緩和時間と拡散係数を評価できる．ただし，計算の途中で等方分布に近いと仮定して，テンソル
３４２

Lecture Note

3. Collision Models for Large Density Scale Plasmas

やベクトルに関する量をスカラー量に低減している．まずが可能である．

式（１４）における速度緩和係数 $は，
$#$
)"$
&，

Y. Sentoku and T. Taguchi

いくつか結果を示す．図５（a）‐（d）は，初期値としてフラットトップな速度分布を与えた粒子群が時間的に緩和し（２５）ていく様子を示したものである．この図で破線は初期値と総エネルギーが等しいマクスウェル分布，すなわちこのエで与えられる．ここで，ネルギーでの終状態である．この終状態は，Rosenbluth

&$# %
$
(
，
$
& #!
''* %
*

ポテンシャルを評価した標的粒子分布でもある．図よりわ（２６）かるように，徐々に終状態分布に漸近して最終的に一致すである．また，拡散係数 " は，る．なお，緩和の様子は２体衝突モデルで計算した結果とほぼ同じであり，時間的変化もよく再現している．

&#*%"# $ $($，''

次に，もう少し平衡状態からずれた初期条件として，２（２７）温度分布，すなわち，速度の速い成分と遅い成分の混じっである．この拡散係数 " は，式（１４）の揺動力の相関関数とた初期分布にした結果を図６に示す．この場合も，マクス

以下のような関係をもつ．

ウェル分布への漸近の様子は２体衝突モデルとほぼ同じで

．
!$
!$
)
"!
#$
)
)
'
##
&
"%
#%
"!)
#%

あることを確認している．このことから，緩和過程は分布（２８）の形状に敏感ではないので，標的粒子の分布を終状態マク
スウェル分布で評価すれば十分である，という Antonsen

実際のプログラムは，各格子における平均速度と平均エネルギーを計算し，これから計算される速度緩和係数 $の仮定が正しいことがわかる．

と拡散係数 " を使って，速度緩和と揺動力を計算し，速度

ただし，この仮定は等方分布の場合には２温度分布のよに修正を加える，という過程を経て次ステップに進む，とうに非平衡性が多少強くても問題ないが，プラズマ中のいう流れになる．
ここで問題がある．前節で詳しく議論したように，このモデルで重要なのは保存量である．基本的に運動量とエネルギーが保存しなければならないが，Langevin 方程式は個々の粒子に独立した修正を加えるので，そのままでは運
動量やエネルギーは保存しない．このとき，あくまでもブラウン運動による平衡状態への緩和さえ合っていればいいと考えれば，厳密に保存する必要はなく，平均量がある誤差範囲で評価できればいいとも言える．しかし，粒子コードの場合，平衡性を評価する局所領域のサイズが小さく，それに含まれる粒子数が少ないので，大きな揺らぎの発生が避けられない．揺らぎの大きさは，電磁場の揺らぎに直接影響を及ぼし，新たな不安定性や人工的な加熱を引き起こす．そこで，緩和を評価するセルごとに修正前と修正後の運動量とエネルギーを計算し，それらが等しくなるように補正を施す．具体的には，各粒子の速度 "(に与える修正

図５

Langevin モデルによる速度緩和の様子１．

図６

Langevin モデルによる速度緩和の様子２．

量を以下の式で評価する．
"("$$!$$"$ %!)!")
$
!)
""(
!"(#(
)"$
&%
)

&

"
"

!)
&$#$
#
)
!
$
("!
．（２９）
%%
''

ここで，!
%は平均が０で２乗平均が１になるような成分を持つ３次元一様乱数である．また，!
! と "は補正パラメー
タであり，! !"(#!

（３０）##"("!"(%#! "(
!$

（３１）

の２式を満足するように選ぶ．式（３０），（３１）は，相対論的公式ではないので，温度が相対論的な場合には使えない．
しかし温度が非相対論的で平均速度が相対論的な場合には，平均速度の系にローレンツ変換をすれば使用すること
３４３

Journal of Plasma and Fusion Research Vol.90, No.6 June 2014

ビーム伝播のような非等方２速度分布ではうまくいかない［１１］．図７はこれを示した図である．
図７（a）は初期分布で，図７（b）はある程度の時間が経過した後の図である．図７（b）には比較のため，２体衝突モ
デルの結果をグレーの細線で示している．この比較からわかるように，簡易 Langevin モデルは２体衝突モデルよりも緩和が速い．ただし，時間が経過すれば，破線のマクス
ウェル分布に 漸 近す る の で，ビ ー ム 伝 播 問 題 で も簡易

図７

Langevin モデルによる速度緩和の問題点．

Langevinモデルが有効であることは確かである．緩和が速い理由は，図７（a）のように破線の標的粒子分布が初期の２速度分布の中間に位置しているためである．標的粒子分布の平均速度と，初期の２速度分布のピークでの速度との差が実際に衝突するピーク間の速度差より小さくなり，速度の３乗に反比例する衝突周波数において無視できない差違となって現れるのである．
ただし，問題になるのは図７のように２個のマクスウェ
ル分布が十分離れている場合，すなわち各ピークの速度拡がりよりも二つの山の速度差が十分大きい場合である．そ
図８

こで，速度拡がりから二つの速度分布を分離することが可

２群簡易 Langevin モデルによる速度緩和の様子．

能であると仮定すれば，粒子を二つの群に分離して，それぞれの群での緩和と，異なる群間の緩和の３種類を計算すない．ただしエネルギーや運動量の保存のための修正にもれば良いという考えが得られる．ここではその詳細を示す

う一度加算計算が必要なので，若干煩雑である．またビーことはできないので，結果だけを示す．

ム相互作用問題のような非等方性の強い分布に対しては，

図８は，初期分布を２群に分けて緩和計算を行った図で同種粒子でも標的粒子の分類が必要になるなどの問題も生ある．図７（b）と比較してわかるように，緩和が遅くなっじる．

ている．この図には，２体衝突モデルも同時に描いている

最終的には，衝突計算がどの程度重要か，プラズマの非が，隠れて見えない．すなわち，ほぼ一致している．この

平衡性，非等方性がどの程度強いか，などに応じて使い分ことから，ビーム分布の場合も，適切な群に分離できれば，ければよいと思われる．

簡易 Langevin モデルが有効であることがわかる．なお，２

３．
４ まとめ

群間の緩和を考えるときには，２群の相対速度方向に平行

今回は，密度が高い領域のシミュレーションに不可欠なな拡散と垂直な拡散の相違も考慮する必要がある．図８の衝突計算を粒子コード中に導入する技法の詳細について述

計算ではこれも含めている．
問題は，レーザープラズマ相互作用で発生する電子ビー

べた．レーザープラズマの解析では，密度変化が大きく，ムの場合には，最初から２群に分離できないことである．

単に高密度領域における衝突を考えるだけではだめで，低

ただ，ここまで見てきたように標的粒子分布の形状には鈍密度領域とのスムーズな接続がポイントである．

感なのであまり精度よく分離する必要はない．経験的に２体衝突モデルは，運動量とエネルギー保存則を満足さは，粒子群をその平均速度の上位と下位に分離するだけでせるのが比較的容易であり，現在もっともよく使われてい十分である．衝突緩和の問題は粒子コードにとっては付加る．ただし，重みの異なる粒子が混在している場合には保価値に過ぎない．さほど精度よく解く必要はなく，あくま

存性を保つために付加的な計算が必要となる．また，粒子でも粒子の運動に補正を与えればよいのである．

をランダムに２個ずつ選択して衝突させる必要があるた

また，群に分けるのが３個以上だと，相互緩和の計算量

め，どの粒子が近傍に存在するかを確認する前処理も必要が増え，群に分離するのも難しくなる．この辺りはトレーで，粒子コードというフレームワークからやや逸脱した感

ドオフで，２体衝突モデルの計算量を凌駕する場合にはエがある．

ネルギーや運動量の保存が保証されている２体衝突モデル

これに対し，Langevin 方程式を用いたモデルは，あくまに軍配が上がる．どちらがよいかは難しいところだが，エでも格子量で速度揺動を計算しようというものなので，粒ネルギーや運動量を自動的に保存するという意味では２体

子コード本来の構造を保っている．ただし，ここで述べた衝突モデルが良い．しかし，粒子の識別が必要なランダム

簡易モデルでは，運動量とエネルギー保存に関する相対論

ペアリングに時間がかかると共に，衝突計算は重心系に移的な拡張が不完全な上に，ビーム伝播のように分布関数が動して散乱させてまた戻すという作業をそれぞれの粒子に２峰性をもつ場合には Rosenbluth ポテンシャルを評価す

ついて計算せねばならず，計算手順は多い．これに対して，

る分布の設定が難しい．ただし，粒子間の直接衝突を計算

簡易 Langevin モデルは温度や流速計算のような加算計算するわけではないので，異なる重みが混在することによるは必要だがそれほど複雑ではなく，計算にかかる労力は少問題は発生しない．このあたりも含めて，今後の発展が待３４４Lecture Note

3. Collision Models for Large Density Scale Plasmas［５］Y. Sentoku and A.J. Kemp, J. Comput. Phys. 227, 6846
(2008).
［６］L. Spitzer Jr., Physics of Fully Ionized Gases (Interscience,
New York, 1956) (Chapter 5).
［７］M.N. Rosenbluth et al., Phys. Rev. 107, 1 (1957).
［８］W.M. Manheimer et al., J. Comput. Phys., 138, 563 (1997).
［９］A. Zhidkov：放電研究 50, 2, 31 (2007).
［１０］T. Taguchi et al., Optics Express 18, 2380 (2010).
［１１］H. Sakagami et al., Laser Part. Beams 30, 243 (2012).

たれるところである．

参考文献［１］R. Shanny et al., Phys. Fluids 10, 1281 (1967).
［２］T. Takizuka and H. Abe, J. Comput. Phys. 25, 205 (1977).
［３］Y. Sentoku et al., J. Phys. Soc. Jpn 67, 4084 (1998).
［４］K. Nanbu and S. Yonemura, J. Comput. Phys. 145, 639(1998).

ひこ

千 徳 靖 彦
Vol.９
０
‐
０
５ p.３
０
４を参照．

!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!

やす

た

ぐち

とし

ひろ

田 口 俊 弘摂南大学理工学部電気電子工学科，教授．
専門はレーザープラズマの理論とシミュ
レーション．現在は，粒子コードや粒子‐流体のハイブリッドコードを使って高密度レーザープラズマの研究をしている．プログラムは，グラフィックルーチンから並列処理言語まで，基本的に何でも自作する主義．１
０年前まではフルマラソンを完走するほどの体育会系であったが，最近は来るべき老後に向けて園芸をやり始めている．

!!!!!!!!!!!!!!!!!!!!!!!!!!

とく!!!!!!!!!!!!!!!!!

せん!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

Y. Sentoku and T. Taguchi

３４５!!!!!!!!!!!!!!!!!!!!!!!!!!