37
hp150123 「京」若手人材育成利用 K Junior Researcher Promotion
磁気リコネクションの 3 次元粒子シミュレーション：
電流層におけるプラズマ乱流
Three-Dimensional Particle-in-Cell Simulation of Magnetic Reconnection:
Plasma Turbulence in the Current Layer
藤本桂三1)
, Richard D. Sydora2) Keizo Fujimoto1), Richard D. Sydora2)
1) 北京航空航天大学 空間与環境学院、2) アルバータ大学物理学科 1) School of Space and Environment, Beihang University, China,
2) Department of Physics, University of Alberta, Canada
要旨 「京」コンピュータを用いて、世界最大規模の 3 次元プラズマ粒子シミュレーションを実施し、 磁気リコネクションに伴って磁気中性線近傍に発生するプラズマ乱流の構造とそれが磁気リコ ネクション過程に与える影響を調べた。その結果、乱流強度やパワースペクトルの冪が計算領域 のサイズによって大きく変わることが見出された。これは、計算領域の拡大にともない、より長 波長モードの不安定波動が励起し、乱流の駆動メカニズムが変わるためである。また、プラズマ 乱流に伴い磁気拡散が増大するため、電流層の幅が拡がる。しかし、従来の理論予測とは異なり、 ある上限を超えて拡がることはできず、電子スケールの幅が維持されることが明らかになった。 キーワード：磁気リコネクション, プラズマ, 粒子シミュレーション, 乱流, 磁気拡散, パワース ペクトル Abstract
Using the K computer, we have performed plasma particle-in-cell simulations of 3D magnetic reconnection with the world’s largest scale, and have investigated plasma turbulence arising in the vicinity of the magnetic neutral line and its impact on the reconnection process. We found that the turbulence intensity and the power of the power spectrum density depend strongly on the computational domain size. This is because the larger domain size yields additional unstable modes with longer wavelength, leading to the change in the driving mechanism of turbulence. In accordance with plasma turbulence, the magnetic diffusion is enhanced, so that the current sheet is broaden. However, it is found that, contrary to the previous theory, the width does not increase beyond an upper limit and remains the electron-scale size.
© 2020 Research Organization for Information Science and Technology All rights reserved. Received: 17 Feburary 2020
Accepted: 6 October 2020 Available online: 16 October 2020
38
Keywords：Magnetic reconnection, Plasma, Particle-in-cell simulation, Turbulence, Magnetic diffusion, PSD 1. 研究の背景と目的 磁気リコネクションは、磁気エネルギーをプラズマの運動エネルギーに効率よく変換すること のできる、重要なプラズマ基礎物理過程の 1 つである。局所的な領域で進行するにもかかわらず、 大域的な磁場の構造を変化させ、大規模なプラズマ対流を引き起こすため、惑星磁気圏や天体周 辺プラズマにおける爆発的なエネルギー解放現象を説明できる、有力な物理過程として注目され ている。また、その過程において非熱的粒子が顕著に生成されるため、高エネルギー宇宙線の発 生源としても期待されている。しかし、クーロン衝突がほとんど起きない宇宙プラズマにおいて、 磁場の繋ぎ変わりがどのように発生するのか（磁気散逸問題）、磁気エネルギーがどのようにし てプラズマ運動エネルギーに変換されるのか（エネルギー変換問題）、そして、どのようにして非 熱的粒子が生成されるのかなど、未解明の問題が数多く残されている。これらの問題を解決する 上で、電子とイオン両方の運動論効果を記述できるプラズマ粒子シミュレーションは、有力な手 段となっている。 近年、計算機性能の急速な向上にともない、3 次元空間における大規模粒子シミュレーション が可能になりつつある。これまで、世界のいくつかのグループにおいて、異なるガイド磁場（初 期電流方向の磁場）強度のもとで計算が実施されてきた[1-3]。これらに共通する結果は、2 次元 シミュレーションでは記述できない強い電磁乱流の発生である。（この結果は定性的には人工衛 星や実験室プラズマの観測事実とも矛盾しない[例えば, 4,5]。）特に、我々のグループでは、最も 基本的な磁場配位である反平行磁場（ガイド磁場無し）の場合について研究を進めている。最近、 「京」コンピュータで実施した大規模シミュレーションの結果から、電流方向の計算領域を従来 よりも拡げた場合、新たに長波長モードの不安定性が励起することが明らかになった[6]。しかし、 新たな波動モードの発生により、磁気中性線近傍の乱流構造やそれが磁気リコネクション過程に 与える影響は十分には明らかになっていない。 本研究では、磁気中性線近傍に焦点を当てて大規模シミュレーションの結果を解析し、プラズ マ乱流構造明らかにし、それが磁気リコネクション過程に与える影響を考察する。 2. 計算モデル
従来のプラズマ電磁粒子（PIC: Particle-in-Cell）モデルに適合細分化格子（AMR: Adaptive Mesh Refinement）を適用した独自の電磁粒子モデル（AMR-PIC コード[7]）を用いた。AMR 法は、計 算格子を与えられた分割指標に従って分割もしくは統合することによって、動的に空間解像度を 調節し、効率的に高解像度計算を実現する手法である。さらに、1 格子あたりの粒子数を調節す るため、超粒子の分割・統合も同時に行っている。分割指標には、各格子点上での Debye 長λDeと
39
初期電流方向の電子流速 Veyが用いられており、∆L > 2λDeもしくは Vey > 2VAが満たされた場合に 格子分割を実施する。ここで、∆Lは格子のサイズ、VAはアルヴェン速度である。
超粒子及び電磁場の時間発展は、従来の PIC コードと同様、Buneman-Boris 法および Yee-Buneman 法（陽解法）を用いた[例えば, 8]。ただし、AMR を実装するため、計算データはツリー 構造を成しており各格子点における物理情報はポインタを介してやり取りされる。 初期の平衡解は Harris タイプの電流層とし、磁場を Bx(z)=B0 tanh(z/δ) 、プラズマ数密度を n(z)=n0 sech2(z/δ) + nb tanh2(z/δ)と与えた。ここで、δは電流層の半値幅である。本研究では、δ=0.5λi、 nb=0.044n0（λiはイオンの慣性長）とした。イオンと電子の質量比は mi/me=100、光速は c/VA=27 とした。イオンと電子の温度比は、T0i/T0e=5、Tbi/Tbe=1、Tbe/T0e=1 と与えた。ここで、T0sと Tbsは、 それぞれ、粒子種 s の初期電流層および背景プラズマにおける温度をあらわす。計算領域のサイ ズは、世界最大規模の Lx×Ly×Lz=82λi×41λi×82λiであり、最大格子サイズを∆LB=0.08λi、最小格 子サイズを∆LD=0.02λiとしたため、最大解像度は 4098×2048×4096 ≈ 4×1010となる。また、最大 粒子数は、イオンと電子の合計で~1012程度である。境界条件は、x と y 方向に周期境界、z 方向 に完全導体壁を与えた。なお、本計算には 80TB 程度の計算メモリが必要であるため、「京」コン ピュータ 8192 ノード（65,536 コア）を占有して実施した。本研究では、比較のために、電流方向 の計算領域を小さくした場合の計算（Lx×Ly×Lz=82λi×10λi×82λi）も実施した。 3. 並列計算の方法と効果（性能） 格子と粒子両方に対して領域分割法による並列化を行なった。従来の PIC コードではポアソン 方程式を用いることが多かったが、その解法（緩和法）には複数回の全体通信が必要であるため、 並列効率の低下を招く要因となっていた。そこで、本研究では、ポアソン方程式を一切使用しな い方程式系を採用し、並列効率を改善した（Fig.1）。また、PIC コードでは粒子の空間分布が時間 とともに変化するため、プロセスあたりの粒子数にばらつきが生じ、負荷バランスが崩れるとい う問題があった。本研究では、プロセスごとの分割領域を動的に変化させ、各プロセスが常に同 等の粒子数を計算するように自動調節する手法（適合ブロック法）を考案した[7]。 Fig. 1：超並列 AMR-PIC コードの性能。各 計算に対して経過時間の逆数がプロットさ れている。縦軸は逐次計算の値で規格化さ れているため加速率を意味する。コード は、それぞれ、電荷保存法と適合ブロック を用いたもの（□印）、電荷保存法と固定ブ ロックを用いたもの（△印）、ポアソン方程 式と固定ブロックを用いたもの（×印）を 示す。 並列化率： 99.8%
40
適合ブロック法の開発により、磁気リコネクションのような複雑な現象のシミュレーションに対 して 99.8%の並列化効率を達成した（Fig.1）。さらに、計算のホットスポットである粒子計算部分 の SIMD 化や Prefetch を実施することにより更なる改善が実現した（Fig.2）。
Fig. 2：さらなる最適化を行なった超並列 AMR-PIC コードの性能。赤線は Fujitsu FX1、青線は 「京」での結果を示す。破線は効率 100%の場合をあらわす。 4. 研究成果 磁気リコネクションは、初期に xz 平面内に小さな磁場擾乱を与え、電流方向（y 方向）に一様 な磁気中性線を発生させることによってスタートさせる。計算開始後しばらくすると、初期の磁 気中性線近傍でテアリング不安定性が局所的に発達し、準定常的な高速磁気リコネクションが始 まる。Fig. 3 に、磁気中性線近傍におけるパワースペクトル密度とリコネクション効率の時間発 展を示す。ここで、パワースペクトル密度は、電場 Eyの単位空間あたりのエネルギースペクトル 密度を示す。また、リコネクション効率は、単位時間に y 軸に沿った単位長さで再結合（リコネ クト）する磁気フラックスの量で定義され、ER≡<-Ey>nline/(VA0B0)で与えられる（<>nlineは平均的な 磁気中性線の位置での y 軸方向の平均値、添え字 nline は磁気中性線を意味する。VA0は上流のア ルヴェン速度）[例えば, 9]。時刻 tωci=13（ωciはイオンの旋回周波数）付近で高速磁気リコネクシ ョンが始まり、これに伴って磁気中性線近傍の波動活動が活発になることがわかる。特に、tωci=15 以降では、電流層から磁気島（フラックスロープ）が間歇的に放出されるため波動活動がより活 発になる[3]。先行研究[6]では、「京」を用いた大規模シミュレーションで電流方向（y 方向）の計 算領域を従来よりも大きくした結果、kyλi≈1.2 の電流層シアー不安定性（CSSI: Current Sheet Shear Instability）[10]に加え、より長波長（kyλi≈0.3）の電子 Kelvin-Helmholtz 不安定性（eKHI: electron Kelvin-Helmholtz Instability）が励起することが明らかになった。
Strong scaling
Effective parallelism 99.9990%
41
Fig. 4 に電場 Eyのパワースペクトル密度（時間平均値）を示す。Ly=41λiとした大規模シミュレ ーション（赤線）では、3 点において明確な冪の変化が見られる。1 つ目は eKHI が励起する kyλi≈0.3、 2 つ目は CSSI の kyλi ≈1.2、そして 3 つ目は kyλi ≈10、すなわち、低域混成ドリフト不安定性（LHDI: Lower Hybrid Drift Instability）の波数である。計算領域が小さな場合（黒線）と比較すると、kyλi <10 において乱流スペクトルが大きく異なっているのがわかる。特に、慣性領域（1.2 < kyλi < 10） における冪が、Ly=10λiの場合は Kolmogorov スケーリング（-5/3）とよく一致するのに対し、Ly=41λi の場合はより急峻になる（スペクトルがソフトになる）ことが明らかになった。 これは、eKHI（kyλi≈0.3）の発生によって、乱流へのエネルギー注入が増加するためであると考 えられる。さらに、この低周波モードの発生によって、乱流の空間当たりのエネルギー密度も大 幅に増加することがわかった。このことから、磁気中性線近傍の乱流過程は、計算領域のサイズ に大きく依存することが示唆される。 Fig. 4：磁気中性線近傍における 電場 Eyの波数空間 kyにおけるパ ワースペクトル密度（計算時間全 体における平均値）。赤線が Ly=41λiの場合で、黒線が Ly=10λi の場合である。黒い破線は、 PSD ∝ (kyλi)α (α = -2.1, -1.6, -4.0)の 関係を示す直線である。 Fig. 3：（a）磁気中性線近傍に おける電場 Eyの波数空間 kyに おけるパワースペクトル密度 (Ey PSD)と（b）リコネクショ ン効率 (ER)の時間発展。
42 電磁乱流により磁気中性線近傍の電流が阻害されるため、磁気散逸・拡散が発生する。衝突系 で一般的に成り立つオーム則を仮定すると、電磁乱流による実効的な異常電気抵抗はηan ≡-<δ(neVe)× δB>y/(<ne><Jy>)と定義される。ここで、任意の物理量 A に対して、A=<A>+δA（<A>は 空間平均値、δA は各点における平均値からの偏差）と表記した。異常磁気拡散の発生に伴い、電 流層の幅が拡がることが期待される。実際、電流層の半値幅δは抵抗長λan≡ηan/(µ0Vin)（µ0は真空 の透磁率、Vinはプラズマのインフロー速度の大きさ）を用いて、以下のように書けることが知ら れている[11]。 𝜹𝜹 𝝀𝝀𝒆𝒆= 𝟏𝟏 𝟐𝟐 𝝀𝝀𝒂𝒂𝒂𝒂 𝝀𝝀𝒆𝒆 + � 𝟏𝟏 𝟒𝟒 � 𝝀𝝀𝒂𝒂𝒂𝒂 𝝀𝝀𝒆𝒆� 𝟐𝟐 + 𝟏𝟏 (𝟏𝟏) ここで、λeは電子の慣性長である。この式から、異常抵抗の増加に従って電流層幅が単調増加す ることがわかる。Fig. 5 に高速磁気リコネクション（ER>0.09）時におけるシミュレーション結果 を示す。図から明確にわかるように、計算領域が大きい場合（赤い正方形）の方が、異常抵抗が 大きくなる。これは、領域サイズとともに乱流強度が増加することと整合的である。一方、理論 曲線と比較すると、異常抵抗が比較的小さい（ηan<0.2）場合は、理論とよく一致するが、異常抵 抗が大きくなると理論曲線からはずれ、δ/λe≈2のほぼ一定値をとることが明らかになった。つま り、どれだけ乱流が強くなろうとも、乱流による磁気散逸・拡散の増大が他の効果によって打ち 消され、ある一定以上には増大しないことがわかった。このことから、強い乱流状態であっても、 電流層の幅は電子スケールよりも大きくなることは無く、乱流による大幅なリコネクション効率 の上昇[例えば, 12]は期待できない。 以上のように、本研究では、電流層に強いプラズマ乱流が発生した場合の、磁気中性線近傍の リコネクション過程が、弱い乱流の場合とは大きく異なることを明らかにした。この強い乱流は、 計算領域を従来よりも拡げたために発生したものである。本シミュレーションを実施するために Fig. 5：磁気中性線近傍における異常抵抗ηanと 電流層の幅δの関係。赤い正方形が Ly=41λiの 場合で、青い正方形が Ly=10λiの場合である。 実線は(1)に基づく理論曲線を示す。
43 は 80TB 程度の計算メモリが必要であり、「京」以外の計算機で実施するのは容易ではない。した がって、本研究成果は、大容量メモリを擁する「京」だからこそ達成することができたと言える。 5. まとめと今後の課題 効率的なエネルギー解放（プラズマ加速）を可能にする磁気リコネクションは、宇宙・天体プ ラズマ現象を解明する上で鍵となる物理過程であるが、その駆動源である磁気散逸機構は十分に は理解されていない。本研究では、大規模な 3 次元プラズマ粒子シミュレーションを実施するこ とにより、磁気中性線近傍に発生する乱流構造とそれが磁気散逸過程に与える影響を調べた。そ の結果、乱流強度やパワースペクトル密度の冪が計算領域のサイズによって大きく変わることが 見出された。これは、計算領域の拡大にともない、より長波長モードの不安定波動が励起し、乱 流の駆動メカニズムが変わるためである。一方、電流層の幅は、乱流強度が非常に大きくなった 場合でも、ある一定値（δ/λe≈2）を上限としてそれ以上増大しないことがわかった。これは、乱流 による磁気散逸の増幅が他の効果によって打ち消されているためである。今後は、計算データの 解析をさらに進めることにより、乱流による磁気散逸を打ち消す効果を同定するとともに、なぜ 大局的な磁気散逸（つまり、リコネクション効率）がある上限値以上に増大できないのかを解明 する。 参考文献
[1] H. Che, J. F. Drake, and Μ. Swisdak (2011), Nature, 474, 184. [2] W. Daughton, et al. (2011), Nature Phys., 7, 539.
[3] K. Fujimoto, and R. Sydora (2012), Phys. Rev. Lett., 109, 265004 [4] M. Zhou, et al. (2009), J. Geophys. Res., 114, A02216.
[5] H. Ji, et al. (2008), Geophys. Res. Lett., 35, L13106. [6] K. Fujimoto (2016), Geophys. Res. Lett., 43, 10557. [7] K. Fujimoto (2011), J. Comput. Phys., 230, 8508.
[8] C. K. Birdsall, and A. B. Langdon (1991), Plasma physics via computer simulation, IOP, London. [9] K. Fujimoto (2009), Phys. Plasmas, 16, 042103.
[10] K. Fujimoto, and R. Sydora (2017), J. Geophys. Res., 122, 5418. [11] V. M. Vasyliunas (1975), Rev. Geophys., 13, 303.
