Particle-in-Cell法

Particle-in-Cell (PIC、セル内粒子) 法とは、特定の問題における偏微分方程式を解く方法の1つである。この方法では、個々の粒子 (または流体要素) が連続な相空間で追跡される。一方で、密度や電流といった分布の積率や電磁場が計算格子上で計算される。粒子の追跡はラグランジュ描像で、積率と電磁場の計算はオイラー描像で記述される。

PIC法は1955年には既に使用されていた^[1]。これは最初のFortranコンパイラが利用可能になるよりも昔の事である。当時の方法は、Buneman（英語版）、Dawson（英語版）、Hockney、Birdsall、Morseらにより、1950年代後半から1960年代初頭にかけてプラズマシミュレーションで人気を博した。プラズマシミュレーションにおけるPIC法では、固定格子上で計算されたセルフコンシステントな電磁場 (または静電場) 内の荷電粒子の軌道を追跡する^[2]。

技術的側面

Buneman、Dawson、Hockney、Birdsall、Morseらにより発明された古典的なPIC法は、多種の問題について直感的であり、かつ簡単に実装できるという長所を持つ。これらの長所があったために、特にプラズマシミュレーションにおける分野でPIC法が成功したと考えられている。古典的なPIC法には、通常、次の手順が含まれる。

運動方程式の積分
電荷と電流の、場の格子への分配
格子上での場の計算
格子から粒子の位置への場の内挿

平均場を介した粒子間相互作用のみを含むモデルは PM (Particle-Mesh) と呼ばれ、直接的な二体相互作用を含むモデルは PP (Particle-Particle) と呼ばれる。また、それらの両方を含むモデルは PP-PM もしくは P³M と呼ばれる。

PIC法は初期の頃から、いわゆる 離散粒子ノイズ による誤差の影響を受けやすい、と認識されてきた^[3]。この誤差は本質的には離散粒子が持つ統計的性質に起因するものである。オイラー法やセミラグランジュ法などの従来の固定格子による手法と比べると、今でも理解が進んでいるとは言い難い。

現代の幾何学的PIC法のアルゴリズムは、非常に多くの理論的枠組みに基づいている。これらのアルゴリズムでは、離散多様体、微分形式の補間、正準または非正準のシンプレクティック数値積分法の手法を使用して、ゲージ不変性、電荷保存則、エネルギー-運動量保存則が保証されると同時に、さらに重要な粒子-場系の無限次元シンプレクティック構造も保証される^[4]^[5]。これらの優れている点は、幾何学的PIC法のアルゴリズムが場の理論の基本的な枠組みに基づいて構築されており、完全形式、つまり物理学の変分原理と直接結びついている事である。

PICプラズマシミュレーション技術の基礎

プラズマの研究では、様々な粒子種（電子、イオン、中性の原子、分子、ダスト粒子など）の系について調査される。PICコードに関連する方程式の組には、ローレンツ力を含む運動方程式と、電場(電界)および磁場(磁界)を解くためのマクスウェル方程式がある。運動方程式を解くコードとマクスウェル方程式を解くコードは分かれており、それぞれ particle mover (または pusher) 、field solver と呼ばれる場合もある。

超粒子

多くの場合、取り扱う系には膨大な数の粒子が含まれており、実行不可能な計算量となってしまう。シミュレーションを効率的に行うために、いわゆる「超粒子 (super particle、または macro-particle)」が使用される。超粒子とは、多数の現実の粒子が集約された、計算上で扱う粒子の事である。例えばプラズマシミュレーションの場合、1つの超粒子は数百万の電子またはイオンに対応する。また、流体シミュレーションの場合、超粒子は渦要素などに対応する。ローレンツ力から受ける加速は電荷質量比のみに依存するため、たとえ超粒子の数をリスケーリングしたとしても、超粒子が現実の粒子と同じ軌跡を辿るようにする事が可能である。

超粒子に対応させる現実の粒子の数 (これを 超粒子の重み という) は、超粒子の運動について十分な統計を収集できるように決める必要がある。系内の異種粒子間 (例えばイオンと中性粒子の間) に大きな密度差がある場合、粒子種ごとに別々の重みを適用する場合もある。

粒子の追跡

超粒子を用いる場合でも、通常シミュレートする超粒子の数は10⁵個以上と非常に多い。運動については各粒子について個別に計算する必要があるため、多くの場合においてPICシミュレーションで最も時間がかかる部分は particle mover コードである。従って、この部分は高い精度と速度をもつ必要があり、様々なスキームの最適化に多大な労力が費やされている。

運動方程式を解く際に使用されるスキームは、陰解法と陽解法の2つに分類される。陰解法 (例えば修正オイラー法) では、同じ時間ステップ内で更新される場の情報から粒子速度を計算するが、陽解法では、前回の時刻における力の情報のみを使用するため、ソルバーは比較的単純で高速に動作する。ただしその代わりに、陽解法では小さい時間ステップ幅が必要となる。PICシミュレーションでは、リープ・フロッグ法がよく使用される。これは、2次の陽解法に分類される^[6]。また、ローレンツ力のうちの電場による加速と磁場による加速を半ステップごとに分離して計算するボリス法もよく使用される^[7]^[8]。

プラズマへ適用させた場合の典型的なリープ・フロッグ法は、次の形式をとる:

{\frac {{\boldsymbol {x}}_{k+1}-{\boldsymbol {x}}_{k}}{\Delta t}}={\boldsymbol {v}}_{k+1/2}

{\frac {{\boldsymbol {v}}_{k+1/2}-{\boldsymbol {v}}_{k-1/2}}{\Delta t}}={\frac {q}{m}}\left({\boldsymbol {E}}_{k}+{\frac {{\boldsymbol {v}}_{k+1/2}+{\boldsymbol {v}}_{k-1/2}}{2}}\times {\boldsymbol {B}}_{k}\right)

ここで、添え字 $k$ は前回の時刻における量である事を示し、 $k+1$ はその次の時刻における量である事を示す(つまり、 $t_{k+1}=t_{k}+\Delta t$ である)。速度は、通常の時刻 $t_{k},t_{k+1},...$ ではなく、それらの中間の時刻 $t_{k-1/2},t_{k+1/2},...$ で計算される。

速度 ( ${\boldsymbol {v}}_{k+1/2}$ 等) はボリス法から求める事ができる。典型的なボリス法は次の形式をとる:

{\boldsymbol {x}}_{k+1}={\boldsymbol {x}}_{k}+{\Delta t}\,{\boldsymbol {v}}_{k+1/2}

{\boldsymbol {v}}_{k+1/2}={\boldsymbol {u}}'+q'{\boldsymbol {E}}_{k}

ここで、

{\boldsymbol {u}}'={\boldsymbol {u}}+({\boldsymbol {u}}+({\boldsymbol {u}}\times {\boldsymbol {h}}))\times {\boldsymbol {s}}

{\boldsymbol {u}}={\boldsymbol {v}}_{k-1/2}+q'{\boldsymbol {E}}_{k}

{\boldsymbol {h}}=q'{\boldsymbol {B}}_{k}

{\boldsymbol {s}}=2{\boldsymbol {h}}/(1+h^{2})

q'=\Delta t\times (q/2m)

である。

ボリス法は長期的に優れた精度を持つため、荷電粒子を追跡する際の事実上の標準となっている。非相対論的なボリス法が長期的に優れた精度を持つ理由としては、シンプレクティック性を持たないにもかかわらず、相空間の体積が保存されるためである事が分かっている。一般にシンプレクティック多様体上のハミルトンフローの性質はボリススキームにも当てはまり、プラズマのマルチスケール（英語版）ダイナミクスの解析には効果的である。また相対論的なボリス法についても、相空間の体積が保存するように修正する事で、交差した電磁場の中で一定速度の解が得られる事が分かっている^[9]。

電場と磁場の計算

マクスウェル方程式 (または、より一般的な偏微分方程式) を解くためによく使用される方法は、次の3つのいずれかに分類される。

有限差分法では、連続領域が点の離散格子に置き換えられ、電場および磁場が計算される。微分は隣接格子点における値の差で近似されるため、偏微分方程式は代数方程式に変換される事になる。

有限要素法では、連続領域が要素の離散格子に分割される。偏微分方程式は固有値、固有ベクトルおよび固有空間として扱われ、最初に、各要素に局所的な基底関数を使用して試行解が計算される。次に、必要な精度に達するまで最適化する事で最終解が得られる。

スペクトル法も、偏微分方程式を固有値問題に変換する。ただし、ここでは基底関数が高次であり、かつ領域全体で定義される。この場合、領域自体は離散化されず、連続領域のまま扱われる。あとは有限要素法と同様に、固有値方程式に基底関数を代入する事で試行解が見つかり、最適化する事で最終解が得られる。

粒子と電磁場の重み付け

「particle-in-cell」という名前は、プラズマのもつ巨視的な物理量 (数密度、電流密度等) が粒子に割り当てられている事に由来する。粒子は連続領域上の任意の位置を取り得るが、一方で巨視的物理量は電磁場と同じように格子点でのみ計算される。そのため巨視的物理量を得るためには、格子点への「粒子の重み付け」を行う必要がある。そこで、1つの粒子はある「形」を持っていると考えて、その「形」が次の形状関数で定められているという仮定を置く:

S({\boldsymbol {x}}-{\boldsymbol {X}})

ここで、 ${\boldsymbol {x}}$ は粒子の位置であり、 ${\boldsymbol {X}}$ は格子点の位置である。形状関数には通常、最も単純で簡単な、一次 (線形) の重み付けスキームが選択される。この手法は、いわゆるセル内クラウド (cloud-in-cell、CIC) 法と呼ばれる。どのようなスキームであれ、形状関数は空間の等方性、電荷の保存、および高次項の精度向上 (収束性) の条件を満たす必要がある^[10]。

電場と磁場は格子点でのみ計算されるため、粒子に作用する力の計算には直接使用できない。そのため、「電場と磁場の重み付け」によって内挿する必要がある:

{\boldsymbol {E}}({\boldsymbol {x}})=\sum _{i}{\boldsymbol {E}}_{i}S({\boldsymbol {X}}_{i}-{\boldsymbol {x}})

ここで、添え字 $i$ は格子点のラベルである。粒子に作用する力がセルフコンシステントに得られるためには、粒子から格子点への数密度および電流密度を割り当てる方法と、格子点から粒子への電磁場の補間方法との間に矛盾があってはならない。なぜなら、いずれの物理量もマクスウェル方程式の中に現れているからである。特に、電磁場の補間スキームでは運動量が保存される必要がある。これは、粒子と電磁場に同じ重み付けスキームを選択し、かつ適切な空間対称性を保証する (すなわち、自己力が無く、作用・反作用の法則を満たす) 事で実現可能である。

衝突

電磁場の計算では自己力を含まない事が必要であるため、セル内では粒子の生み出す電磁場が粒子から遠ざかるにつれて小さくなる必要がある。すると、セル内での粒子間力が過小評価されてしまうが、これは荷電粒子間のクーロン衝突（英語版）を加える事でうまくバランスを取る事ができる。大きな系で全ての組の相互作用を考慮しようとすると計算負荷が高くなりすぎる事から、代わりにいくつかのモンテカルロ法が開発された。その中でも広く使用されるものは2体衝突モデルである^[11]。2体衝突モデルでは、同一セル内にある粒子の中から無作為にペアが抽出され、ペアの粒子同士が衝突する。

実際のプラズマでは、荷電粒子-中性粒子間の弾性衝突、電子-中性粒子間の電離衝突などを含む非弾性衝突から化学反応に至るまで、クーロン衝突以外の衝突も重要な要素となる場合があり、各衝突は区別して扱われる必要がある。荷電粒子-中性粒子間の衝突を処理する衝突モデルの殆どは、全ての粒子の衝突確率を直接計算するDSMC法か、荷電粒子種ごとの最大衝突確率で処理するヌル衝突法を使用する^[12]^[13]。

精度と安定性の条件

他のシミュレーション手法と同様に、PIC法でも時間ステップと格子サイズの幅を適切に選択する必要がある。これらは、興味のある現象が時間と長さのスケールで適切に解かれるようにするだけでなく、コードの処理速度と精度にも影響する。

陽的時間積分スキーム (例えば、広く普及しているリープフロッグ法等) を使用した静電プラズマシミュレーションに対しては、解の安定性を確保するために、時刻ステップ幅 $\Delta t$ および格子サイズ幅 $\Delta x$ について、次の2つの重要な条件を満たす必要がある。

\Delta x<3.4\,\lambda _{D}

\Delta t\leq 2\,\omega _{pe}^{-1}

これらの条件は、非磁化の1次元プラズマ調和振動を考えると導かれる。後者の条件は厳密に満たされる必要があり、シミュレーションを実施する上では、エネルギーを保存させる目的でより厳しい制約が要求される。そこで係数2の部分は1桁小さい数値に置き換えられ、一般的には $\Delta t\leq 0.1\,\omega _{pe}^{-1}$ が使用される^[14]。このように、プラズマの自然な時間スケールは逆プラズマ周波数 $\omega _{pe}^{-1}$ によって、長さスケールはデバイ長 $\lambda _{D}$ によって与えられる。

また、陽的な電磁プラズマシミュレーションでは、時間ステップ幅はCFL条件も満たす必要がある。

\Delta t<\Delta x/c

ここで、 $\Delta x\sim \lambda _{D}$ であり、 $c$ は光速である。

応用

プラズマ物理学の分野では、PICシミュレーションは、レーザー-プラズマ相互作用、オーロラを発生させる電離層内の電子加速とイオン加熱、磁気流体力学、磁気リコネクション、トカマク中のイオン温度勾配不安定性やその他の微視的不安定性、真空アーク（英語版）、およびダストプラズマの研究で役立っている。

純粋なPIC法だけでなく、流体モデルとのハイブリッドモデルが使用されることもある。ハイブリッドモデルでは、いくつかの粒子種の運動の処理にはPIC法が使用され、他の粒子種は流体モデルでシミュレートされる。流体モデルで扱う粒子種は、特定の速度分布(通常はマクスウェル-ボルツマン分布)に従うと仮定される。

また、PICシミュレーションはプラズマ物理学以外の固体力学および流体力学の分野でも利用されている^[15]^[16]。

PICアプリケーション

商用製品

製品名	Web	国内販売元 (開発元)
LSP	^[17]	株式会社エーイーティー (Alliant Techsystems Inc.)
MAGIC	^[18]	株式会社エーイーティー (Alliant Techsystems Inc.)
Particle-PLUS	^[19]	株式会社ウェーブフロント (株式会社ウェーブフロント)
VSim	^[20]	株式会社エーイーティー (Tech-X Corporation)

オープンソースまたは研究用のコード

ソースコード名	Repo	ライセンス	研究論文
ALaDyn	^[21]	GPLv3+
EPOCH	^[22]	GPL	doi:10.1088/0741-3335/57/11/113001
FBPIC	^[23]	三条項BSD-LBNL	doi:10.1016/j.cpc.2016.02.007
iPic3D	^[24]	APL 2.0	doi:10.1016/j.matcom.2009.08.038
open-pic	^[25]	MIT License	doi:10.3847/1538-4357/aa6d13
piccante	^[26]	GPLv3+
PICLas	^[27]	GPLv3+	doi:10.1016/j.crme.2014.07.005 doi:10.1063/1.5097638
PIConGPU	^[28]	GPLv3+	doi:10.1145/2503210.2504564
PicUp3D		GPLv2	doi:10.1109/TPS.2006.883402
SHARP			doi:10.3847/1538-4357/aa6d13
Smilei	^[29]	CeCILL-B	doi:10.1016/j.cpc.2017.09.024
UPIC	^[30]
VPIC	^[31]	三条項BSD	doi:10.1063/1.2840133
Warp	^[32]	三条項BSD-LBNL	doi:10.1063/1.860024
WarpX	^[33]	三条項BSD-LBNL	doi:10.1016/j.nima.2018.01.035
XOOPIC	^[34]		doi:10.3847/1538-4357/aa6d13
ZPIC	^[35]	AGPLv3+

教科書

Charles K. Birdsall and A. Bruce Langdon, 'Plasma Physics via Computer Simulation', McGraw-Hill (1985), ISBN 0-07-005371-5
Roger W. Hockney and James W. Eastwood, 'Computer Simulation Using Particles', CRC Press (1988), ISBN 0-85274-392-0

外部リンク

テキサス大学講義: Computational Physics - Particle-in-cell codes

脚注

^ F.H. Harlow (1955). A Machine Calculation Method for Hydrodynamic Problems. Los Alamos Scientific Laboratory report LAMS-1956.
^ Dawson, J.M. (1983). “Particle simulation of plasmas”. Reviews of Modern Physics 55 (2): 403–447. Bibcode: 1983RvMP...55..403D. doi:10.1103/RevModPhys.55.403.
^ Hideo Okuda (1972). “Nonphysical noises and instabilities in plasma simulation due to a spatial grid”. Journal of Computational Physics 10 (3): 475–486. Bibcode: 1972JCoPh..10..475O. doi:10.1016/0021-9991(72)90048-4.
^ Qin, H.; Liu, J.; Xiao, J.; et al. (2016). "Canonical symplectic particle-in-cell method for long-term large-scale simulations of the Vlasov-Maxwell system". Nuclear Fusion. 56 (1): 014001. arXiv:1503.08334. Bibcode:2016NucFu..56a4001Q. doi:10.1088/0029-5515/56/1/014001。
^ Xiao, J.; Qin, H.; Liu, J.; et al. (2015). "Explicit high-order non-canonical symplectic particle-in-cell algorithms for Vlasov-Maxwell systems". Physics of Plasmas. 22 (11): 12504. arXiv:1510.06972. Bibcode:2015PhPl...22k2504X. doi:10.1063/1.4935904。
^ Birdsall, Charles K.; A. Bruce Langdon (1985). Plasma Physics via Computer Simulation. McGraw-Hill. ISBN 0-07-005371-5
^ Boris, J.P. (November 1970). "Relativistic plasma simulation-optimization of a hybrid code". Proceedings of the 4th Conference on Numerical Simulation of Plasmas. Naval Res. Lab., Washington, D.C. pp. 3–67.
^ Qin, H.; et al. (2013). "Why is Boris algorithm so good?" (PDF). Physics of Plasmas. 20 (5): 084503. Bibcode:2013PhPl...20h4503Q. doi:10.1063/1.4818428。
^ Higuera, Adam V.; John R. Cary (2017). “Structure-preserving second-order integration of relativistic charged particle trajectories in electromagnetic fields”. Physics of Plasmas 24 (5): 052104. Bibcode: 2004JCoPh.196..448N. doi:10.1016/j.jcp.2003.11.004.
^ Tskhakaya, David (2008). “Chapter 6: The Particle-in-Cell Method”. In Fehske, Holger; Schneider, Ralf; Weiße, Alexander. Computational Many-Particle Physics. Lecture Notes in Physics 739. 739. Springer, Berlin Heidelberg. doi:10.1007/978-3-540-74686-7. ISBN 978-3-540-74685-0. http://cds.cern.ch/record/1105877
^ Takizuka, Tomonor; Abe, Hirotada (1977). “A binary collision model for plasma simulation with a particle code”. Journal of Computational Physics 25 (3): 205–219. Bibcode: 1977JCoPh..25..205T. doi:10.1016/0021-9991(77)90099-7.
^ Birdsall, C.K. (1991). “Particle-in-cell charged-particle simulations, plus Monte Carlo collisions with neutral atoms, PIC-MCC”. IEEE Transactions on Plasma Science 19 (2): 65–85. Bibcode: 1991ITPS...19...65B. doi:10.1109/27.106800. ISSN 0093-3813.
^ Vahedi, V.; Surendra, M. (1995). “A Monte Carlo collision model for the particle-in-cell method: applications to argon and oxygen discharges”. Computer Physics Communications 87 (1–2): 179–198. Bibcode: 1995CoPhC..87..179V. doi:10.1016/0010-4655(94)00171-W. ISSN 0010-4655. https://zenodo.org/record/1253854.
^ Tskhakaya, D.; Matyash, K.; Schneider, R.; Taccogna, F. (2007). “The Particle-In-Cell Method”. Contributions to Plasma Physics 47 (8–9): 563–594. Bibcode: 2007CoPP...47..563T. doi:10.1002/ctpp.200710072.
^ Liu, G.R.; M.B. Liu (2003). Smoothed Particle Hydrodynamics: A Meshfree Particle Method. World Scientific. ISBN 981-238-456-1
^ Byrne, F. N.; Ellison, M. A.; Reid, J. H. (1964). “The particle-in-cell computing method for fluid dynamics”. Methods Comput. Phys. 3 (3): 319–343. Bibcode: 1964SSRv....3..319B. doi:10.1007/BF00230516.
^ “ATK - LSP”. Northropgrumman.com. 2020年11月12日閲覧。
^ “ATK - MAGIC”. Northropgrumman.com. 2020年11月12日閲覧。
^ “WaveFront - Particle-PLUS”. Wavefront.co.jp. 2020年11月12日閲覧。
^ “Tech-X - VSim”. Txcorp.com. 2020年11月12日閲覧。
^ “ALaDyn”. Github.com. 2020年11月12日閲覧。
^ “EPOCH sign in”. Cfsa-pmw.warwick.ac.uk. 1 December 2017閲覧。
^ “FBPIC”. GitHub.com. 2020年11月12日閲覧。
^ “iPic3D”. GitHub.com. 2020年11月12日閲覧。
^ “XOOPIC”. GitHub.com. 2020年11月12日閲覧。
^ “piccante”. GitHub.com. 2020年11月12日閲覧。
^ “PICLas”. GitHub.com. 2020年11月12日閲覧。
^ “PIConGPU”. GitHub.com. 2020年11月12日閲覧。
^ “Smilei”. GitHub.com. 2020年11月12日閲覧。
^ “UPIC”. GitHub.com. 2020年11月12日閲覧。
^ “VPIC”. GitHub.com. 2020年11月12日閲覧。
^ “Warp”. Bitbucket.org. 2020年11月12日閲覧。
^ “WarpX”. GitHub.com. 2020年11月12日閲覧。
^ “XOOPIC”. Langmuir.eecs.berkeley.edu. 2020年11月12日閲覧。
^ “ZPIC”. GitHub.com. 2020年11月12日閲覧。