無衝突衝撃波のダイナミクス
加藤 恒彦∗ （国立天文台）

1

無衝突衝撃波と宇宙の高エネルギー現象
無衝突衝撃波とは、粒子間の衝突がほとんど起きない無衝突プラズマ中に発生する衝撃波であり、衝撃波における散逸過程は粒子間の衝突ではなく、そこで作られる電場や磁場を介して行われる。宇宙空間でよく見られる高温で希薄なプラズマは無衝突プラズマであり、そこで発生する衝撃波は無衝突衝撃波であると考えられる。
無衝突衝撃波は宇宙における高エネルギー現象、特に非熱的な高エネルギー粒子の生成に対して非常に深い関係がある。まず、衝撃波の散逸が粒子間の衝突ではなく電磁場によるものであるため、電場や磁場をうまく利用できた一部の粒子が選択的に高エネルギーにまで加速される可能性がある。また、このような加速機構自体ではそれほど高エネルギーにまで加速できない場合でも、ある程度の「中間」のエネルギーにまで加速された粒子が、フェルミ加速などのさらに高エネルギーにまで加速できる加速機構の「種粒子」になる可能性がある。また、衝撃波の下流での「熱的」な粒子の分布も Maxwell 分布になるとは限らず、電子・陽子プラズマなどではそれぞれの下流温度も一般に異なる。このことは、衝撃波の Jump condition に影響し、間接的にフェルミ加速の効率にも影響を与えることになる。
さらに、近年ガンマ線バーストやそのアフターグロウの衝撃波に対して考えられているように [1]、衝撃波面付近で生成された磁場が、衝撃波の下流である程度の時間残る可能性もある。このような磁場は、高エネルギー粒子のシンクロトロン放射に対して重要になり、また、フェルミ加速における散乱体として働いてその加速効率に影響を与える可能性もある。特に、相対論的な衝撃波では、多くの場合、背景磁場の向きが衝撃波面法線に対して垂直になってしまうので1 、そのままではフェルミ加速は働かないが、衝撃波下流で背景磁場と同程度以上の強さの乱れた磁場が存在すれば、それにより粒子が散乱されフェルミ加速が働く可能性がある。ガンマ線バースト、AGN ジェット、パルサー風などの多くの相対論的な衝撃波では、power-law的な高エネルギー電子の存在が観測により示唆されており、もしこれらの電子の加速機構がフェルミ加速であるとすると、このような磁場の揺らぎが衝撃波面付近で生成されているのかもしれない。
以上のような点を詳しく知るためには、無衝突衝撃波自体の物理を知ることが不可欠である。しかし、その物理過程は電磁場とプラズマを構成する粒子の運動論的な振る舞いが複雑に関係する非線形現象であるため、もっとも有力な研究手段は、コンピュータによるシミュレーションである。以下では、電子・陽電子プラズマ中の相対論的な無衝突衝撃波の２次元シミュレーションの結果を紹介する。特に、(1) 衝撃波面付近で磁場が作られるのか、(2) 衝撃波の散逸はどのように行われるのか、(3) 衝撃波付近で発生する電磁場に付随した粒子加速機構があるか、に着目する。シミュレーションは、背景磁場と衝撃波の構造との関係を見るため、初期に背景磁場が無い場合とある場合の２通りについて行う。

∗ Email: katoutn@cc.nao.ac.jp

1 上流の静止系で上流磁場が任意の向きを向いていても、電磁場のローレンツ変換を考えると衝撃波静止系では垂直成分がローレンツ因子倍になるためである。

2

無衝突衝撃波のシミュレーション

2.1

電磁粒子シミュレーション法

以下で行うシミュレーションは電磁粒子シミュレーションと呼ばれるもので、電磁場と個々の粒子の時間発展を consistent に解くものである [2, 3, 4]。このタイプのシミュレーションでは Particle in Cell と呼ばれる方法を用いる。これは、電磁場はグリッド上で Maxwell 方程式を解いて求め、粒子については個々の粒子の運動を相対論的な運動方程式を解いて求めるものである（図 1 参照）。この方法では粒子の運動量空間の分布もきちんと取り入れられるため、さまざまな運動論的な現象を扱うことができる。一方で、空間的には電子の skin depth や Debye length が、時間的には電子のプラズマ振動数やサイクロトロン振動数の逆数が基本的スケールとなるので、その適用は主にプラズマのミクロスケールの現象に限られる。

図 1: シミュレーション法
以下ではシミュレーションの単位として次の量を用いている。
−1/2[s]−1/2[cm]−1• 時間： τ0 = ωpe∼ 1.8 × 10−5 ne0−1• 長さ： l0 = cωpe

∼ 5.3 × 105 ne0• 電磁場： E∗ = B∗ = c(4πne0 me )1/21/2∼ 3.2 × 10−3 ne0[esu] or [gauss]

ここで ωpe ≡ (4πne0 e2 /me )1/2 は電子プラズマ振動数であり、ne0 は考えている系の特徴的な電子数密度である。最右辺は cgs-Gauss 単位系での値である。また、粒子の速度は主に４元速度 u を用いて表す。通常の速度 v と対応するローレンツ因子 γ = (1 − v 2 /c2 )−1/2 を用いると u = γv/c である。

2.2

シミュレーションの設定

シミュレーションの設定は、背景磁場を除いて２つのシミュレーションで共通である（図 2 参照）。
組成は電子・陽電子プラズマである。計算領域は２次元で、シミュレーション平面を x − y 平面とする。

物理サイズは Lx × Ly = 480 × 60、グリッド数は Nx × Ny = 4096 × 512 である。粒子は電子、陽電子それぞれについて約 3 億 2000 万粒子を使用する。

衝撃波は、粒子を反射する壁にプラズマをぶつけることにより発生させる。初期状態ではプラズマ粒子を２つの反射壁の間に一様に配置し、各プラズマ粒子には平均として右向き（+x 方向）に ux = 2.0 の初速度を与えておく。計算をスタートさせると、壁で反射された粒子と流入してくる粒子との間でプラズマの不安定性が発生し、衝撃波が形成されてゆく。なお、シミュレーションの座標系は衝撃波の下流系に対応し、下流静止系で衝撃波が右から左へ伝播していく様子を見ることになる。

図 2: シミュレーションの設定

2.3

背景磁場が無い場合の衝撃波

背景磁場が無い場合の衝撃波のシミュレーションでは、初期の電場・磁場は系全体にわたってゼロにしている。
計算終了時（t = 233）の電子数密度を図 3 に示す（下流系で見た上流プラズマの数密度を１としている）。
左から上流のプラズマが右へ流れてきて、x = 350 − 400 の間の非一様な領域を経て x > 400 のほぼ一様な下流領域に至ることがわかる。

図 3: 背景磁場が無い場合の衝撃波の電子数密度
図 4 は、y 方向に平均した数密度、x 方向の平均速度、電場・磁場の大きさのプロファイルである。数密度と平均速度は、いずれも x = 350 − 400 の領域で遷移していることがわかる。マクロに見れば、この領域が「衝撃波面」に相当する。そして重要な点は、この遷移領域で強い磁場ができていることである。この磁場のエネルギー密度は上流の粒子のバルクの運動エネルギーの約 8% に達し、後に見るように衝撃波の散逸は主にこの磁場が担っている。さらに、衝撃波の遷移領域の上流にも比較的強い電場と磁場が存在していることがわかる。
図 5 に、磁場（Bz ）と電流密度（Jx ）を示す。衝撃波面付近で作られた強い磁場の変動の方向（波数ベクトルの方向）は主に y 方向であることがわかる。この磁場の構造を作っているのは、x 方向に流れる多数の電流のフィラメントであり、この図からも磁場の構造と電流密度のパターンが対応していることがわかる。

図 4: y 方向に平均した電子の数密度（左）、x 方向の平均速度（中央)、電場・磁場の大きさ（右）のプロファイル

図 5: 磁場 Bz （上）と電流密度 Jx （下）衝撃波面付近で流れの方向の電流フィラメントが多数でき、それが磁場を作り出すという状況は、ガンマ線バーストの衝撃波などで提案された Weibel 不安定性による磁場生成モデル [1] と consistent である。Wiebel不安定性とは、プラズマ粒子の速度分布が非等方的な場合に発生する無衝突プラズマ中の不安定性であり[5, 6]、不安定性の進行に伴って数多くの電流フィラメントができ、それにより磁場が作られる（図 6 の左を参照）。電流フィラメントは、その中の電流が Alfvén 電流 [7] に達するまで成長を続け、サチュレーション時の磁場は、初期の非等方性が大きい場合には粒子の運動エネルギーの数パーセントのエネルギーに達する[8]。２つのプラズマの成分が counter streaming 的（２流体的）になっている場合は速度分布の非等方性の極端な場合であり、このような場合にも Weibel 型の不安定性が発生して電流フィラメントおよび磁場が作られることが粒子シミュレーションにより示されている [9, 10] （図 6 の右を参照）。このような Weibel 型の不安定性でできる磁場の波数ベクトルは、速度分散の大きい方向に垂直、あるいは counter streaming plasmaの場合は流れの方向に垂直である。今回のシミュレーションの設定では、波数ベクトルの方向は y 方向になり、これはシミュレーション結果の磁場構造と consistent である。以上により、衝撃波面で発生している電流および磁場は、Weibel 型の不安定性により生成されたと考えてよいと思われる。なお、今回の２次元シミュレーションでは系の対称性から磁場が Bz 成分しか作られず、したがって衝撃波面において粒子の速度の z 成分は熱化されない点に注意する必要がある。

図 6: Weibel 不安定性
図 7 は、電子の位置 x と４元速度の x 成分 ux の位相空間プロットである。上流からの電子が衝撃波遷移領域で一気に熱化されている様子がわかる。特に、衝撃波面付近で一部の粒子が下流から上流へ向かって逆流している点が重要である。この逆流粒子が上流から流入してくる粒子との間で counter streaming 的な状態を作り、それが Weibel 不安定性を発生させ、磁場が作られると考えることができる。またここでは示さないが、同時に静電ビーム不安定性2 も発生しており、それによる静電場も衝撃波面の少し上流側で見られる。さらに、これらの逆流粒子は、衝撃波面のさらに上流側で磁場の擾乱を作る。この磁場に上流からの電子と陽電子がぶつかることで荷電分離・電流生成が起き、図 4 に見られる比較的強い磁場と電場が上流にできると考えられる。
電子のエネルギー分布を図 8 に示す。左の図は、下流領域（x = 430 − 442）の電子のエネルギーのヒストグラムである。参考のため、３次元の場合の相対論的 Maxwell 分布と２次元の場合の相対論的 Maxwell 分布も示しているが、ヒストグラムはどちらにも良くは一致しない。すでに述べたように２次元シミュレーションの対称性から粒子は z 方向には熱化されないので、３次元の Maxwell 分布に一致しないのは当然であるが3 、２次元の Maxwell 分布にも一致しないのは、衝撃波の散逸過程においては磁場による等方化が主であり、エネルギーの緩和が十分に起きないためであると考えられる。
2 同様に counter streaming 的な状況で発生する不安定性で、正味の電荷密度が発生して静電場ができる。流れの方向に波数ベクトルを持つ縦波モードである。
3 ３次元で同様のシミュレーションをすれば当然各方向に熱化される。３次元の電子・陽電子の衝撃波のシミュレーションは最近
Haugbølle や Hededal らによって行われたが [11]、下流の粒子の分布関数についての詳細な解析はまだされていない。

図 7: 電子の x − ux 位相空間プロット
図 8 の右の図は系全体の電子のエネルギー分布である。power-law 的な高エネルギー粒子の生成は特に見られない。

図 8: 電子のエネルギー分布：下流領域の電子の分布（左）と系全体の電子の分布（右）2.4背景磁場がある場合の衝撃波

背景磁場がある場合の衝撃波のシミュレーションでは、背景磁場およびそれに伴う電場を初期に与える。
それ以外の点では背景磁場が無い場合と同様である。
初期の磁場は y 方向（流れの方向、したがって衝撃波面法線の方向に垂直）に与え、その強さは σ パラメータにより与える。σ パラメータとは、衝撃波系で見た上流の磁場のエネルギーフラックスと粒子のエネルギーフラックスの比として定義されるもので、Crab nebula のモデルでよく用いられる [12, 13]。しかし、衝撃波の伝播速度はシミュレーションの結果わかるものであるため、ここでは、下流系で見たエネルギーフラックスの比として σ パラメータを定義する。σ パラメータの値については、Kennel & Coroniti による Crabnebula の MHD モデル [12, 13] で最もよく観測を説明できるとされている値 σ = 0.003 を使用する。
なお、電場は上流静止系でゼロになっていると考えるのが自然である。シミュレーション系（下流系）では上流プラズマが動いているため、motional な電場が −z 方向に発生する。初期の磁場を B0z とすると、この電場は E0z = −(Vx /c)B0z で与えられる（Vx は初期のプラズマのバルク速度である）。初期に系全体にこの電場を設定する。

図 9 は計算終了時（t = 303）の電子数密度である4 。σ = 0.003 程度の磁場であっても、衝撃波の構造は背景磁場が無い場合とは大きく異なることがわかる。特に、衝撃波遷移領域の上流側において、密度の揺らぎなどの構造は特に無く、遷移領域の幅も背景磁場が無い場合に比べて小さくなっている。

図 9: 背景磁場がある場合の衝撃波の電子数密度
図 10 は、y 方向に平均した数密度、x 方向の平均速度、y 方向の磁場のプロファイルである。x = 360 − 380の領域で数密度と平均速度が遷移しているのがわかる。磁場 By は、これよりも少し広い領域にわたって激しく変動し、最終的には一定値に収束する。この系での衝撃波の伝播速度は −0.32c であり、衝撃波系で見たときの下流磁場と上流磁場の強さの比を見積もると約 3.0 となる。これは MHD 近似 [12] が与える値 2.96とよく一致する。

図 10: y 方向に平均した電子の数密度（左）、x 方向の平均速度（中央)、y 方向の磁場（右）のプロファイル
衝撃波面付近の磁場の構造を図 11 に示す。磁場は、単に By が圧縮されるだけでなく、Bx や Bz 成分も衝撃波遷移領域で上流磁場の 10 倍程度の強さに達することがわかる。Bx 、Bz の磁場の波数ベクトルの方向は y 方向であり、背景磁場が無い場合と同様、電流フィラメントの生成によりできた構造と思われる。Byは、遷移領域で激しく変動し、その向きが反転する領域も見られる。なお、Bz で見られる格子状のものは、おそらく数値的な原因による電磁波と思われる。
図 12 の左の図は、y 方向に平均した磁場で規格化した磁場の揺らぎ（平均からの差）の大きさを示す。これを見ると、衝撃波遷移領域の下流においても、平均磁場の半分ぐらいの強さの磁場の揺らぎが残っていることがわかる。これは右の図からわかるように、衝撃波下流で z 方向の電流が残るためである。この電流がどれくらいのタイムスケール生き残ることができるのかについてはまだ良くわからないが、もし長期的に生き残って磁場の揺らぎを維持するのであれば、Crab nebula の観測から示唆されている衝撃波下流での乱れた磁場 [14] の起源になっている可能性もある。
4 x > 420 に見られる構造は境界の影響と思われる。時間が経っても、この領域の構造はほとんど変化しない。

図 11: 磁場の構造：Bx （左上）、By （左下)、Bz （右上）

図 12: 磁場の揺らぎq∆Bx2 + ∆By2 （左）と z 方向の電流密度 (右）なお、図 12 の左の図では、衝撃波の上流側にもある程度強い磁場の細かな揺らぎが存在していることがわかる。これは、Gallant らの１次元シミュレーションにより報告された、衝撃波面における synchrotronmaser 不安定性に起因する coherent な電磁波放射と思われる [15]。電磁波は衝撃波面から上流へ向かって伝播し、precursor emission と呼ばれる。
電子のエネルギー分布を図 13 に示す。左の図は衝撃波下流領域の電子のエネルギー分布である。赤線がシミュレーション結果のヒストグラム、青線が３次元の Maxwell 分布である。両者は完全には一致しないものの、大体においては一致しているといえる。
（なお、今回のシミュレーションでは最初に背景磁場 By がある
ため、粒子は z 方向にも速度を持つことができ、結果として衝撃波面で全ての方向に熱化が起きる。）
図 13 の右の図は、系全体の電子のエネルギー分布である。高エネルギー側において、わずかではあるがpower-law 的な分布が見られる。ただし、その power-law index は 9.4 と非常に steep である。

図 13: 電子のエネルギー分布：下流領域の電子（左）と系全体の電子（右）の分布
すでに述べたように、垂直衝撃波ではフェルミ加速が働くことは難しく、また、仮に働く条件がそろっていたとしても、フェルミ加速は、粒子が磁場の乱れにより散乱されてランダムウォーク的に動く過程で衝撃波面を何度も往復してエネルギーを得ていくメカニズムなので、今回のシミュレーションの空間サイズおよび計算時間の範囲内では十分に加速過程を追えていないと考えられる。一方、このような垂直な背景磁場がある衝撃波では、背景磁場に伴う motional な電場を利用した加速が起こりうる。ドリフト加速やサーフィン加速 [16] などの加速機構でも、基本的にはこの motional な電場の方向に粒子が動くことでエネルギーを得ていく。今の場合、上流の背景磁場に対応して −z 方向に電場が生じている。したがって、もしこのような種類の加速が主に働いているのであれば、ローレンツ因子の増分 ∆γ と電子の z 方向の変位 ∆z との間には次の関係がある。

∆γme c2 = −eEz ∆z
実際には Ez の値は衝撃波遷移領域付近で激しく変動するので、z と γ のプロットが１つの直線上に乗ることはないが、少なくとも両者の間に相関があることは推測できる。図 14 がシミュレーション結果の z − γ プロットである。予想されたとおり、γ が大きい（エネルギーが高い）電子は、z の値も大きい傾向があることがわかる。この加速は、ドリフト加速またはサーフィン加速と同様なものかもしれないが、その詳細は今後の解析で明らかになると思われる。

参考文献[1] M. V. Medvedev and A. Loeb, ApJ 526, 697 (1999).
[2] C.K. Birdsall and A.B. Langdon, Plasma Physics via Computer Simulation, (McGraw-Hill, NewYork,1985).

図 14: 電子の z − γ プロット。なお、初期の電子のバルクエネルギーに対応する平均的ローレンツ因子は√
γ = 5 であり、初期の z の値は 0 である。

[3] 田中基彦、西川恭治, 「高温プラズマの物理学」, （丸善、1991）.
[4] M. Hoshino, PhD thesis, Univ. Tokyo (1986) .
[5] E. S. Weibel, Phys. Rev. Lett. 2, 83 (1959).
[6] B. D. Fried, Phys. Fluids 2, 337 (1959).
[7] H. Alfvén, Phys. Rev. 55, 425 (1939).
[8] T. N. Kato, Phys. Plasmas 12, 080705 (2005).
[9] Y. Kazimura, J. I. Sakai, T. Neubert, and S. V. Bulanov, ApJL 498, L183 (1998).
[10] L. O. Silva, R. A. Fonseca, J. W. Tonge, J. M. Dawson, W. B. Mori, and M. V. Medvedev, ApJL 596,L121, (2003).
[11] T. Haugbølle, PhD thesis, Niels Bohr Institute (2005) .
[12] C. F. Kennel and F. V. Coroniti, ApJ 283, 694 (1984).
[13] C. F. Kennel and F. V. Coroniti, ApJ 283, 710 (1984).
[14] S. Shibata, H. Tomatsuri, M. Shimanuki, K. Saito and K. Mori, MNRAS 346, 841 (2003).
[15] Y. A. Gallant, M. Hoshino, A. B. Langdon, J. Arons and C. E. Max, ApJ 391, 73 (1992).
[16] M. Hoshino and N. Shimada, ApJ 572, 880 (2002).