論 文
θピンチの粒子シミュレーション
―運動論的効果―
坂井一雄 竹内智 松本道男
（平成2年8月31日受理）
Particle Simulational Study of Theta Pinch
-Kinetic
Effects-KazuoSAKAI SatoshiTAKEUCHI MichioMATSUMOTO
Abstract Theta−pinch processes are studied by using a 2 Ei 2・dimensional fully electromagnetic parti− cle code in the cases with the compressing magnetic field and static magnetic field both parallel and antiparallel． It is shown that kinetic effects such as E×B drift and the electric・field− gradient effect are responsible for the compression and heating of plasma in the early phases in both parallel and antiparallel cases． The compression of plasma is accompanied with the appearance of the shock wave． In the case where the compressing field and the static field are antiparallel， the rapid acceleration of electrons are observed in the very early stage． The acceleration is related to the trapping of electrons by the magnetic neutral sheet． §1 序 θピンチによるプラズマの圧縮加熱システムは，安 定性において良好な性質を示し，比較的簡単に優れた 閉じこめパラメタを実現している1・2）。理論的には，対 象とするプラズマがいわゆる高βプラズマで，しかも 経時的変化を本質的に伴うため，その振る舞いの理論 解析はもっぱら磁気流体的手法によってきた3）。一方， 実験的には，無衝突衝撃波の形成や電子の異常加熱が 報告され4・5），磁気流体的解析の限界が明らかになり， 運動論的取扱いが必要になってきた。しかしながら， 経時変化する不均一高βプラズマの運動論的解析は 容易でなく，電子加速の機構はいまだ明らかにされて いない。 このようなプラズマの研究には粒子シミュレーショ ＊電子情報工学科，Department of Electrical Engineering＆ Computer Science ンが有効であると考えられる。われわれの研究の目的 は，強いθピンチ的電磁場とプラズマの相互作用を粒 子シミュレーションにより追跡し，主として，電子加 速の機構を明らかにすることにある。ここで紹介する シミュレーションはまだ試験的な段階にあるが，すで に大変興味ある結果が得られているので，それを報告 したい。特に，われわれが強い波動と荷電粒子との非 線形相互作用の基本過程の一つであることを示してき た，Vp×B加速6−11）との関連についても議論する。
Vp×B加速は静磁場Bに垂直に位相速度Vpで伝
播する電磁波動によって荷電粒子が捕捉され，波動と 共に速度VpでBを横切って運動すると，ローレソツ カqVp×B／cによって，波面に沿って加速される現象 である。静電波の場合6−8），荷電粒子は静電ポテンシャ ルによって捕捉され得ることは容易に理解されよう。 この場合は既に原理実験によって，実際にこの機構に より電子が加速されることが確認されている9）。われわれは，さらに波動が横波の場合にも10川，その波動の 位相速度が真空中の光速cより小さければ，同様の加 速が起こることを示した。横波によるVp×B加速の 本質的特徴は，静磁場により磁気中性面の位置が波動 のみのときの位置からずれることと磁気中性面により 粒子が捕捉されることにある。波動のみのときは磁気 中性面は波動磁場が零になる位置に形成されるが，そ こでは波動電場も零になるので，荷電粒子は捕捉はさ れるけれど，加速されない。しかし，適当な大きさの 静磁場を予め掛けて置けぽ，磁気中性面の位置は波動 電場が零の位置からずれる。したがって，Vpとほぼ等 しい速度で動く荷電粒子はこの中性面に捕捉され，中 性面に沿って存在する横電場によって，磁場に曲げら れずに中性面に沿って加速されることになる。ところ で，θピンチにおいては，磁気中性面が形成されるよう な配位が存在するから，その場合にも，この加速機構 が働く可能性があるものと予想されるのである。 §2 シミュレーション・モデル シミュレーションは標準的な2次元電磁粒子コード を使用して実施された。粒子シミュレーションは，原 理的には，ヴラソフ方程式を生成する変分汎関数をい わゆる超粒子によって構成された試行関数の変分に対 して停留させる解を求めることと等価である。ここで， 超粒子というのは位相空間の中で局所的にのみ値を持 つ分布関数である。実際には，ある時刻tにおける電磁 場中で通常の運動方程式に従う超粒子の運動を微小時 間△tだけ追跡し，新しい時刻t＋△tにおける位置お よび速度を求め，その電荷および電流を空間に張られ た格子点に分配し，電荷密度分布および電流密度分布 を求める。次に，それらが生成する，時刻t＋△tにお ける電磁場をマックスウェル方程式を解いて決定す る。この手続きを繰り返すことによって系の運動論的 振る舞いを追跡することができる12）。 さて，われわれのコードにおいては，電磁場による ローレンツカを受けて運動する粒子の運動方程式は非 相対論的なものを用いた。この運動方程式は数値的に は蛙飛びの方法によって解かれた。また，電磁場の決 定には，複素場13）に対するマックスウェル方程式を フーリエ変換した式が用いられた。この式は粒子の運 動方程式と酷似した形になるので，運動方程式の解法 と同じアルゴリズムが使える。数値的には，空間的に は高速フーリエ変換が，時間的には同じく蛙飛びの方 法が使用された。超粒子の形状因子は幅が格子点の数 倍のガウス分布とした。この形状因子は，短波長成分 をカットすることにより，衝突効果が必要以上に強調 Z ↑
一一一＞X
O
9 ： 昌x
86
名 IBa
periodic
←z
ω 匡 鱈工PLASMA
o
←国parabolic
国 工 レりperiodic
図1 シュミレーション・システム 系の左半分を示す．右半分は中心について点対称である． されるのを防ぐ働きがある。 シミュレーション・システムは図1に系の左半分が 示されている。まず，y軸に平行で，一様な静磁場属 の中に，x座標について放物線的で， z座標について一 様な初期密度分布を持つプラズマがある。プラズマの 初期速度分布はマックスウェル分布で，その初期温度 は一様で，イオン温度T，と電子温度Teは等しいとし た。また，プラズマを構成するイオンと電子の質量比 は人工的な値5を用いた。これは試行的シミュレー ションとして，θピソチの全行程を比較的短時間に計 算させるために採用したものである。他に有効な理論 解析の手法がない現状では，この程度のシミュレー ションでもθピンチの定性的理解の第一歩としては 十分有用だと考えられる。 さて，θピンチ装置ではピンチコイルが設置される が，シミュレーションではプラズマの表面近傍でz軸 に平行に流れるシート電流 ノ＝A ｛exp（−at）−exp（一 bt）｝ （1） によって，θピンチ的変動電磁場をプラズマに印加す る。ここで，tは時間である。Aの符号の正負はそれぞ れ順磁場配位と逆磁場配位に対応する。ここで順（逆） 磁場配位というのはプラズマ表面近傍で静磁場と同じ （逆）向きの変動磁場を生成する電流の向きをいう。 電磁場に対する境界条件はx端で完全導体壁，2端で 周期的境界条件を採用した。 また，初期における種々のパラメタは，電子熱速度 vth。／c＝0．46，静磁場の大きさ疏／4πneλD＝O．56，ωp。／ ωce＝3．88， aλD／c −4．34×10−4，およびb＝5aとした。 ここで，eは電気素量， nは粒子数密度，λDはデバイ長， ωp，は電子プラズマ振動数，またωceは電子サイクロト ロン振動数である。このときアルフェーン速度はVA＝ O． 56 v，h，となる。これらはすべて中心部における初期の 一52一値である。系の大きさは4c＝192ArD，1，＝48λ，である。超 粒子の数は電子とイオンそれぞれ16384個である。電子 熱速度として実際の値よりかなり大きな値が設定され ているが，これはデバイ長やラーモア半径を大きくす ることによって運動論的効果を際だたせるためであ る。このように現段階では，粒子シミュレーションの 手法は写実的ではなく，むしろ印象派（点描画）的で
0
th
t＝90
ある。その役割は実際を完全に模倣するのではなく， 有効な直観的描像を獲得することにあるといえる。 §3 シミュレーションの結果 シミュレーションの結果をまず順磁場配位につい て，ついで逆磁場配位について以下に記述する。3．1 順磁場配位
シート電流の式においてAを正に取ると，変動電磁 場は図2の下部のような空間構造になり，プラズマ中 に磁気中性面は形成されない。また，同図上部に点描 されている電子の位相空間分布，x−v、分布，が示すよう に，電子分布の空間的圧縮と表面近傍での電子の速度 分布の広がりが生じているが，後に示す逆磁場配位に おいて起きるような負の2方向への異常な電子加速 は見られない。この図は時刻t＝90におけるものである が，以後プラズマ全体の圧縮と加熱がが徐々に進行し ていく。なお，時間は特徴的時間λD／cによって規格化 されている。プラズマ表面の位置の時間変化を図3に0
図2 順磁場配位における変動電磁場の空間変化と電子x−・Vz の分布． 口は変動磁場By，●は縦電場Ex，○は横電場島を表 す．横軸κの1目盛りは規格化された値で30，縦軸の190
60
民30
0 0 60 120 180 240 t 図3 順磁場配位における表面と衝撃面の位置の時間変化． ●はプラズマ表面，＋は衝撃面に対応する．直線部の傾き からそれらの移動速度として Vsurface＝027c および VSh，，k・O．40c を得る．t＝0
t＝30
t＝60
t＝90
t＝120
t＝150
t＝180
t＝210
図4 順磁場配位における電子の位相空間分布の時間変化 空間的には系の左半分を示した．左欄はx−Vx分布，右欄 はx−Vz分布である．横軸xの1目盛りは図2と同じであ る． 示した（図4も参照せよ）。これから，表面はほぼ一定 速度vsurface・0．27cで移動することが分かる。また，こ れはイオンの位相空間分布においてより明瞭なのだ が，衝撃面がプラズマ中を伝播していくのが観測され る。この伝播速度は当然表面の移動速度より速く， v、h。，k＝O．40cである（図3）。 どうしてこのようになるかを理解するのはさほど難 しくない。図2に示されているように，変動電磁場は プラズマの内部に入るにつれて急激に減少するが，プ ラズマ表面近傍では，静磁場も変動磁場も正なので全 磁場Bも当然正である。また，電場は，昆が負で，盈 が正である。従って，まず，負のEzによって電子およ びイオンはxの正方向にE×Bドリフトをし，その結 果プラズマは圧縮されるものと解釈される。E，もByも 空間的に激しく変化するので，一様モデルに基づく結 果を厳密には適用できないけれど，1つの目安として その結果を用いることは許されるだろう。そこで，い わゆるドリフト速度を推算してみると，プラズマ表面 でおよそVd＝（E，／By）c −O．27cとなり，表面の移動速 度と良い一致を示す。 次に，図4により位相空間分布をより詳細に検討し てみると，初期においてX−Vz分布の両端が斜めに削い だようになっていることが見て取れる。これは有限 ラーモア半径効果で説明できる。すなわち，時計回り に電子は回転するので，プラズマの左端からラーモア 半径程度までの間には上行する電子しか存在しえない のである。実際削がれた部分の幅はその場における ラーモア半径にほぼ等しいことが分かった。このよう な描像を描くと，右行する電子と左行する電子の数は 等しいはずだから，x−・Vx分布にはこのよう斜削部は存 在しない。しかし，以下に示すように左端では右行す る電子の方が左行する電子より速いのでVx分布は上 方にシフトすることになる。実際図4はそうなってい る。 時間が経過すると表面の速度分布の幅が広がってく る。すなわち，電子の加熱が生じている。この機構は 次のようなものだと考えられる。図2から分かるよう に電場の主成分はEzであり，表面近傍でそれは負で， 勾配dEz／duは正である。また，磁場はyの正方向にか かっている。磁場は実際には空間的に変化しているけ れど，副次的要素で本質的な部分を隠さないために， 一様静磁場であると仮定しよう。このような電磁場の 中で原点を中心にして円運動する電子を考える，図5 のように電子は上行するときカーθ互（−rL）で加速され， およそラーモア半径rLだけ走り，円の上方を右行する とき最も速くなる。一方，下行するときはカーeE。（rL）を 受けてrLだけ走り，減速され，円の下方を左行すると き最も遅くなる。したがって，1回転の間に電子にな される仕事△Wはおよそ △w・erL｛E。（rL）一阜（−r、）｝・2・rL・」巖 （2）
となる。これは正であるので，電子は加熱されること 一54一図5 Ez勾配による加熱機構 時計回りに回転する電子は電場五2により円軌道の左半分 で正の仕事，右半分で負の仕事をされる．上図の場合1回 転当たりの仕事は正になる。 になる。これから単位時間当たりの温度上昇率は
ナ誓舌陽 （3）
この結果は電荷の正負によらないので，イナンにも適 用できる。これがθピンチにおけるプラズマ加熱の主 たる機構であると考えられる。3．2 逆磁場配位
今度は，Aを負に取り，シート電流の向きを逆転さ せると，変動磁場Byはプラズマ表面近傍でづ方向を 向き，Aの絶対値を十分大きく取れぽ，全磁場の値が 零になる面，すなわち，磁気中性面が生じる。このた め逆磁場配位では順磁場配位では考えられないような ことが起こり，現象は複雑になる。図6の上部には電 子の位相空間分布の左半分がt＝60，120および180につ いて点描されている。下部には変動電磁場の空間構造 が示されている。なお，現象は系の中心に関して点対 称なのでこれ以降の記述は左半分について行う。磁気 中性面は一対が磁場の極小点の両側に出現するが，横 波によるVp×B加速理論によれぽ1°），左側の磁気中 性面は，そこではE。が正なので，荷電粒子を捕捉できo
0
0
tε180 ；， 壁ら是 灘‘識纏｛
耀懸騨
；、、 5♪ ilY 0 一B。 o 一B。O
一Bo 図6 逆磁場配位における変動電磁場の空間変化と電子のx−Vzの分布 口は変動磁場疏●は縦電場E。○は横電場Ezを表す．変動磁場を静磁場B。だけバイアスすると全磁場が でる．補捉中性面の位置が矢印で示されている．目盛りは図2と同じ．0．5 N
国
二国
0．00
一H＋一一H＋一十＋＋＋͡＋一什一＋＋＋一＋＋＋＋＋＋＋＋＋＋t ●●●●●●● ・°・る・・9．s 8●●●060000
50
100
．．・…°°°°°・°° O OOOO Ot l50
60
30
閣 三． R ●●°°．●・ × §Ov
200
図7 捕捉中性面の位置とその位置における電場の時間変化 X，rapは捕捉中性面の位置座標である．縦の波線は磁気中性面形成の時刻を示す．磁気中性面形成以前におい ては磁場の極小点の位置座標Xmi。とそこにおける電場の値を描いてある。＋はXmmまたはX，．ap，●はE。， ○はEzを表す． ず，反発的に作用する。これを反発中性面と呼ぽう。 一方，右側は荷電粒子の捕捉に有効である。この捕捉 に有効な磁気中性面は矢印で指示されている。これを 捕捉中性面と呼ぽう。 初期（t＝30）においては電子分布の左端の膨張（−x 方向への移動）が観測される。これは，静磁場と逆向 きの印加磁場により全磁場の絶対値が小さくなるため に，ラーモア半径が伸長することによる。半径が引き 伸ばされた電子の中にはシート電流の裏側まで達する ものが現れる。磁気中性面が出現するまでは，全磁場 が極小になる位置はほぼ一定で，x＝30位の所にある。 その点におけるEzは正で，値は0．2程度に達する（図 7）。したがって，初期分布の表面付近の電子はこのE， により一z方向に加速されると共に，小さくなった磁場 によって時計回りに大半径で回転させられる。このた めに，x−Vz分布の左端はV字型になる。また，一部は シート電流の裏側に入って，壁にまで達する。壁では 計算上粒子は完全反射される。壁付近では磁場そのも のも，その勾配も大きいので，これらの粒子は壁に沿っ てドリフトする。 磁気中性面が形成された（t；60）後に，捕捉中性面 の位置で一2方向への急速な電子加速が生起している のが観測される。このような一z方向への異常な電子加 速は順磁場配位では見られなかったもので，磁気中性 面の存在と密接に関連している。 図7に示されているように，捕捉中性面は生成後t＝ 90頃までは，平均速度0．3c程度で右側へ移動するが， その後しぼらく停滞し，t・−150頃から再び右側への移 動を再開する。同図に示されているように，t＝90頃ま では捕捉中性面上では瓦が正なので，これによって電 子は一z方向に磁場に曲げられずに急速に加速された のだと解釈することができる。電子は加速されるとま すます強くその中性面に捕捉されるので，被加速粒子 の位相空間分布は図6でも分かるように横幅が狭い楕 円形の塊になる。このような初期には捕捉中性面上に 負のExも働き，これは被捕捉電子が捕捉中性面ととも にXの正の方向に移動するとき正の仕事をし，瓦とと もに電子加速に寄与する。 捕捉中性面が停滞中にはかなり強い正の島が中性 面上に生じている。しかし，被捕捉粒子も当然停滞す るので，、Exは仕事をせず，被捕捉電子の加減速には影 響を与えない。従って，この停滞期にはわずかに存在 する正のEzによって，被捕捉電子の加速がゆっくりと 進行する。この停滞期にはプラズマの圧縮は進行しな い。これは変動電磁場の印加とともにすぐ圧縮が始ま る順磁場配位との大きな違いである。t・150を過ぎる と，プラズマは圧縮され始め，中性面も再びxの正の 方向に移動し始める。この時期には中性面上に正の瓦 一56一と正のExが存在し，これらは被捕捉粒子に符号が逆の 仕事をするため，被捕捉電子の加速は飽和してしまう。 さて，反発中性面の右側にあった電子は右へとさら に反発され，捕捉中性面に収束していく。一方，反発 中性面の左側にあった電子はさらに左へと反発され， 壁に押し付けられる。ところで，時間経過と共に反発 中性面付近の星は大きくなり，かつdEz／dUも正なの で，順磁場の所で述べたように，壁に押し付けられた 電子はこのEz勾配によって加熱される。壁付近の縦に 非常に細長い分布はこのようにして生じる。 壁に押し付けられた集団と捕捉中性面に捕捉された 集団を除いた残りのバルク電子の分布は順磁場の場合 とよく似た形状をしていることに注意せよ。圧縮が進 行すると共にバルク電子は順磁場の場合と同様にEa 勾配によって加熱される。 §4 結 語 以上の粒子シミュレーションの結果から，実際のθ ピンチ過程においても次の運動論的効果が働くのでは ないかという示唆が得られる。まず順・逆両磁場配位 どちらでも有効な効果として （1）E×Bドリフトによる圧縮効果 （2）横電場勾配による加熱効果 があり，逆磁場配位においてのみ有効な効果として （3）ラーモア半径伸長による膨張効果 （4）磁気中性面による粒子捕捉と加速効果 がある。（3），（4）の効果はθピンチ過程の初期に働き， （4）の効果は，いわゆる横波によるVp×B加速であり， 予想どうり急速な異常電子加速をもたらす。また，E× ．Bドリフトによる急激な圧縮は衝撃面の形成に導く。 ところで，今回のシミュレーショソは試行段階なの で，それなりの限界がある。まず，われわれのモデル ではプラズマの形状は2次元矩形であった。しかし， 実際は円筒状であろう。この差は衝撃面の伝播速度や 表面の移動速度など定量的結論には影響を与えるだろ うが，上述の効果の有効性そのものを否定することは ないだろう。つぎに，人工的な質量比5，電子熱速度 O．46cという値が用いられ，比較的短時間で圧縮加熱 が進行しているので，種々の不安定性や衝突効果など に比して，運動論的効果が強調され過ぎている嫌いは ある。より大きな質量比とより小さな電子熱速度に対 するシミュレーションを行い，いろいろなパラメタに よって結果がどのように変わるか吟味する必要はあろ う。そうするためには長時間の計算が必要になるので， 長大な計算を許す環境の成熟を待たねばならない。今 後の課題としたい。
