J. Plasma Fusion Res. Vol.93, No.6 (2017)２
８５‐２８６

講座

実践低温プラズマシミュレーション
Practical Low-Temperature Plasma Simulation

１．はじめに
1. Introduction
小田昭紀
ODA Akinori
千葉工業大学（原稿受付：２
０
１
７年５月１
１日）

電離度が低く，かつ熱的に非平衡（プラズマ中の正イオ

マシミュレーションを実施し，得られた結果を可視化する

ンおよび原料ガスの温度に比べて電子温度が圧倒的に高までの一連の流れを連続的に行えるようになってきておい）な性質をもつ低温プラズマを利用した技術は，古くはり，低温プラズマシミュレーションを実施する敷居が昔よ

蛍光灯，その後，半導体製造プロセス，産業用光源やプラりも格段に低くなってきている．

ズマディスプレイパネル，オゾン生成など，現在では医しかしながら，この低温プラズマシミュレーションを，療・農業水産分野へとその応用分野を広げながら利用され自分で行うにせよ，商用ソフトウェアを使用して行うにせており，我々の生活を支える基盤技術となっている．

よ，そのシミュレーションがどういった前提（仮定）や方

この低温プラズマ技術を開発するにあたり，多くの研究法で行われており，何をどの程度知ることができるのか

者により理論や実験などを通じて精力的な研究開発がなさを，ある程度正しく理解しておかねば，誤った条件で得られてきたことは今さら言うまでもないところであるが，それたシミュレーション結果を（悪気無く）正しい結果と解の際，この低温プラズマのもつ特性を把握する手段とし釈して使用されてしまうことが強く懸念される．

て，様々な計測手法に基づくプラズマ診断および計算機シ

本講座では，本学会誌で出版された数多くの低温プラズ

ミュレーションが挙げられる．特に，低温プラズマの計算

マシミュレーションに関する小特集・講座・解説記事［６‐機シミュレーション，すなわち低温プラズマシミュレー１０］とは趣が異なり，低温プラズマシミュレーションに興

ションは，上述のプラズマ技術の中でも特に半導体プラズ

味のある方（学生・院生・技術者・研究者）を対象に，核マプロセス技術の進展に伴い１
９８０年代から急速に発展し，融合や高温プラズマのモデル化との相違を明確にしなが現在では低温プラズマの特性を解明する手段として欠かす

ら，シミュレーションを実施することの意義，何がどの程ことのできないものになっている．それに際しては，低温度わかるのか，実験とどの程度一致するのかなどの様々なプラズマを適切に数値モデル化し，それに基づく各種支配観点から低温プラズマシミュレーションの現状をできるだ

方程式を初期条件ならびに境界条件の下で数値的に計算をけ平易に順序立てて紹介する．本講座を通してご一読いた実施する．その際，昔は低温プラズマシミュレーションをだくことで，低温プラズマシミュレーションの興味が広が実施する者自身がプログ ラ ミ ン グ 言 語（例 え ば，FOR-り，更なる理解を促すきっかけとなることを願ってやまな

TRAN や C など）を使用してシミュレーション用プログラい．

ムをコーディングしていたのに対し，最近では自分でコーディングすることなく，低温プラズマシミュレーションを

本講座の構成は，以下の通りである．

行うための商用ソフトウェア［１‐５］が数多く出てきてお

第１章

はじめに

小田昭紀（千葉工大）

り，あたかも PC 上で図を描画するかのようにモデル化を

第２章

低温プラズマシミュレーションの基礎行い，その上で初期条件や境界条件を設定した後にプラズ

２．
１ 低温プラズマとはChiba Institute of Technology, Narashino, CHIBA 275-0016, Japanauthor’s e-mail: akinori.oda@it-chiba.ac.jp

２８５!2017 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.93, No.6 June 2017２．
２ プラズマと固体・液体表面（界面）における反応（Particle- In-Cell/Monte Carlo collision method），中性粒子

２．
３ 低温プラズマシミュレーションの重要性の挙動にDSMC法（Direct Simulation Monte Carlo method）

池田

圭（㈱アテナシス）を適用した各種材料プロセス用プラズマシミュレーション

２．
４ 粒子モデルによるシミュレーション手法の事例を紹介する．３．
２節では，流体モデルによるプロセス

２．
５ 流体モデルによるシミュレーション手法

プラズマシミュレーションと題し，代表的な低温プラズマ中舘
第３章

博，松永史彦（ペガサスソフトウェア㈱）応用の１つである半導体製造プロセスで多く利用されてい低温プラズマのシミュレーション事例

るプラズマエッチングおよびプラズマ CVD 装置を適切に

３．
１ 粒子モデルによるプロセスプラズマシミュレーションモデル化した低温プラズマシミュレーション事例について

田中正明（ペガサスソフトウェア㈱）紹介する．その際，得られた結果からシミュレーションに

３．
２ 流体モデルによるプロセスプラズマシミュレーション

より何がどの程度わかるか，実験との比較を含め紹介す池田圭（㈱アテナシス）

る．３．
３節では，流体モデルによる大気圧プラズマシミュ

３．
３ 流体モデルによる大気圧プラズマシミュレーション

レーションと題し，従来の材料プロセスのみならず医療・小田昭紀（千葉工大）
第４章

バイオ応用へと更なる応用展開がはかられている大気圧低

おわりに温プラズマのシミュレーションの一例として，誘電体バリ
田中正明（ペガサスソフトウェア㈱）ア放電型低周波駆動された大気圧プラズマのモデル化，実験との比較を通じたモデルの妥当性検証，その上での本プ

第１章では，本講座の趣旨を説明し，次章以降の章立て

ラズマシミュレーションから得られた成果を順に解説すを紹介する．

る．

第２章では，本講座で対象とする低温プラズマとは何か

第４章では，本講座を総括するとともに，低温プラズマ

を，核融合や高温プラズマの特性（プラズマ中の粒子の種シミュレーションをとりまく課題等に関して述べる形で本類，電離度や電子，イオン・ガス温度の関係など）との比講座を締めくくる．

較を行いつつ説明する．その上で，低温プラズマシミュ
レーションとは何か，そして，なぜシミュレーションを行

参考文献

う必要があるのかを説明する．その後，低温プラズマシ［１］PEGASUS; http://www.psinc.co.jp/products/［２］CFD-ACE+; http://athenasys.co.jp/main/product/cfdace/ cfd-ace.html［３］ACE+ Suite; https://www.esi-group.com/jp/softwaresolutions / virtual-environment / cfd-multiphysics / acesuite［４］COMSOL Multiphysics; http://www.kesco.co.jp/comsol/［５］VizGlow; http://www.wavefront.co.jp/CAE/VizGlow/［６］真壁利明，南部健一：プラズマ・核融合学会誌 77, 1007
(2001).
［７］浜口智志：プラズマ・核融合学会誌 80, 110 (2004).
［８］斧 高一：プラズマ・核融合学会誌 80, 909 (2004).
［９］杤久保文嘉：プラズマ・核融合学会誌 91, 307 (2015).
［１０］杤久保文嘉：プラズマ・核融合学会誌 92, 680 (2016).

ミュレーションのモデル化手法である粒子モデルおよび流体モデルに関し解説する．
第３章では，各種用途で利用される低温プラズマのシ
ミュレーション事例の紹介を行う．３．
１節では，粒子モデル
によるマグネトロンスパッタプラズマシミュレーションと題し，プラズマ材料プロセスにおける低温プラズマシミュ
レーションには，流体モデル，ハイブリッドモデル，そして粒子モデルが使用されるが，装置条件や運転条件に応じて各モデルの特長を考慮してモデルを選択することが重要である．本節では，粒子モデルによる低温プラズマシミュ
レー シ ョ ン の 例 と し て，荷 電 粒 子 の 挙 動 に PIC/MC 法!!!!!!!!!!!!!!!!!!!!

お

だ

あき

のり

小 田 昭 紀２
００
１年北海道大学大学院工学研究科電子情報工学専攻博士後期課程修了．博士（工学）．現在，千葉工業大学工学部教授．低圧
および大気圧条件下の低温プラズマ（弱電離非平衡プラズマ）の基礎特性解明に関する研究に従事．これまで本学会誌にて数回各種記事を執筆しておりましたが，２
０
１
５年７月から本学会誌の編集委員を拝命しております．当学会誌読者が今回の講座記事を一読いただくことによって，低温プラズマそれ自身，そして低温プラズマシミュ
レーションに関して少しでも興味が湧いてもらえれば非常に嬉しく思います．

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

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

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

２８６

J. Plasma Fusion Res. Vol.93, No.6 (2017)２
８７‐２９１講座 実践低温プラズマシミュレーション

２．低温プラズマシミュレーションの基礎
2. Fundamentals of Low-Temperature Plasma Simulation
池田圭１），中 舘博２），松 永 史 彦２）

IKEDA Kei１）, NAKADATE Hiroshi２）and MATSUNAGA Fumihiko２）
１）

株式会社アテナシス，２）ペガサスソフトウェア株式会社（原稿受付：２
０
１
７年５月１
１日）

低温プラズマは弱電離でかつ熱的に非平衡であり，それゆえに利用価値が高く，様々な産業に幅広く応用されている．しかしながら，実用化されている技術が数多くある一方で現象の理解は必ずしも進んでいない．本章
では，はじめに温度に着目した上で低温プラズマの特徴を説明する．次に，低温プラズマがどのような産業で応用されているかを鳥瞰し，その理解のためにはシミュレーションが不可欠な理由を述べる．更に，低温プラズマ
シミュレーションの手法としての粒子モデルおよび流体モデルに関して，その特徴と基礎的事項を説明する．
Keywords:
low-temperature plasma, non-equilibrium plasma, particle model, PIC/MC, DSMC, fluid model

２．
１ 低温プラズマとは温度は一桁以上異なっていることが多く，非平衡性は維持

我々の身近にあるプラズマの多くは，電離している粒子

されていることが本学会誌の過去の小特集でも取り上げらの数が相対的に少ない弱電離プラズマである．中性粒子がれている［３］．低温プラズマと高温プラズマについては，何室温であると仮定して，弱電離プラズマ中のイオンや電子

度以上は高温プラズマといった厳密な区別はないが，身近の温度はどうなっているであろうか．電子は質量が小さくな低温プラズマは熱的に非平衡であり，その多くは重粒子

電界からの加速を受けやすい一方で，一回の衝突で失うエの温度が気体と同程度の温度から熱電離が大きく影響を及

ネルギーは重粒子同士の衝突の場合よりも非常に小さい．

ぼさない 5,000 K 程度までである（アークの条件を除くと，半導体製造プロセスで多く利用されるような低圧下では，多くの低温プラズマの電離度は 0.1％ にも満たない）．大気

粒子同士の衝突周波数も低いために，プラズマ中で加速さ圧程度の高圧力で重粒子の温度がおよそ 10,000 K まで高くれた電子の温度の方が重粒子の温度よりも高い状態［１］となると熱プラズマと呼ばれるようになり，より熱平衡に近なり，熱的に非平衡なプラズマ，すなわち低温プラズマがづくが，依然として非平衡性が議論されている［４］．更に温容易に生成される．

度が上昇して電離度も高くなると高温プラズマと呼ばれ，衝突周波数は圧力に比例して高くなることから，図１に核融合プラズマはその典型である．２．
２節では身近な低温示すように，大気圧近傍もしくはそれ以上の圧力では，電プラズマに焦点を絞り，どのような分野で応用されている子の温度も重粒子の温度に近づくと期待される［２］．しかのか，また２．
３節では，なぜシミュレーションが必要とされしながら，実際には大気圧でも電子温度と周囲の重粒子のるのかを順に述べる．

２．
２ プラズマと固体・液体表面（界面）における反応
プラズマを扱う装置のシミュレーションを検討する際には，プラズマ領域のモデル化は勿論であるが，周囲の壁との相互作用の考慮も重要である．反応場をごく一部の限られた領域で制御できる場合でも，ガスや電力などを供給するシステムが直接的・間接的に影響を与える．プラズマを応用する目的は，大きく分けて以下の２つに分けられる．
一つは，プラズマを利用した燃焼［５］のように，主にプラズ
マ生成部を利用する場合であり，もう一つはプラズマを固
図１

体や液体［６］と接触させて，その界面との相互作用を利用

電子温度および重粒子の温度の圧力依存性（文献［２］
より
引用）
．

する場合である．後者の代表的な例の一つである半導体の

1)ATHENASYS Co.,Ltd., 2)PEGASUS Software Inc.

authors’ e-mail: dk@athenasys.co.jp, h-nakadate@psinc.co.jp, matsunaga@psinc.co.jp

２８７!2017 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.93, No.6 June 2017

製造過程には，薄膜を堆積するプラズマ CVD（化学気相成た領域で進展し，光電離が重要なプロセスの場合には，そ長）や微細加工を行うエッチングの工程が多数あり，LSIの反応を測定する手段もない．更に，表面反応に関しては（大規模集積回路）や LCD（液晶ディスプレイ）などのデバ

汎用のモデル構築も容易ではない．例えば，文献等で見つイス作製に必須の技術となっている［７］．デバイス応用以けられる金属表面の二次電子放出係数は，一般に研磨され外でも幅広い応用が進んでおり，例えば DLC（ダイヤモン

た状態で測定された値であるが，多くの低温プラズマにおドライクカーボン）［８］等のコーティング技術は，耐摩耗性いて，その表面状態は履歴に依存しやすく，表面の堆積物を必要とする研磨剤や切削工具，自動車部品の他，PETもしくはエッチングされた状態を考慮する必要がある．

ボトルのバリア膜［９］などに応用されている．また，親水性

以上のように，低温プラズマの反応は極めて複雑である向上を目的とした表面改質の例として，大気圧プラズマにため，現象の理解にはシミュレーションによる可視化が不よるチタンインプラントの表面処理も知られている［１０］．

可欠である．以降，具体的な低温プラズマシミュレーショ

異種材料の接着性を向上させる例として，有機膜とステンンの方法を粒子モデル，流体モデルの順に説明する．

レスの密着性向上の目的にも応用が期待されている［１１］．

２．
４ 粒子モデルによるシミュレーション手法

更に，自立式鉄塔としては世界最高となった東京スカイツ
リーの完成には溶接技術が不可欠であった［１２］．その他に

低温プラズマシミュレーションでは，荷電粒子（電子おも，生成物を直接利用する代わりに，プラズマアクチュ

よびイオン）だけでなく，中性粒子（原料ガス，励起原子・エータ［１３］を利用した気流制御技術の研究も進められてお分子，ラジカル分子）を取り扱う必要がある．荷電粒子にり，再生可能なエネルギーの一つとして風車のブレードへ

ついてはの応用も期待されている［１４］．このように，小さなスケー（PIC/MC 法）［２２，２３］，中性粒子については Direct Simula-Particle-in-Cell/Monte Carlo Collision method

ルから大きなスケールまで，実に幅広い分野で利用・研

tion Monte Carlo method（DSMC 法）［２４，２５］という粒子モ

究・応用がなされている．

デルのシミュレーション手法がある．どちらも短い時間ス
テップ t で粒子の運動を追跡していく手法である．

また，近年では医療応用に対する研究も進められており，２００９年 に は ISPM（International Society for Plasma

両者の粒子の移動や境界条件の扱いはよく類似していMedicine）［１５］も設立され，外傷や手術の際にプラズマをる．計算コードのソースレベルの話をすれば，PIC/MC

照射して止血を始めとする創傷医療やがん治療［１６，１７］，法と DSMC 法の粒子境界条件，粒子の移動，統計処理など

心筋梗塞の緩和［１８］治療などへの研究も進められている．

のルーチンは共通化して完全に同一のルーチンを用いるこ

この他にも，種子の発芽制御や成長促進，滅菌・殺菌や鮮とができる．
超粒子（super particle や simulated particle などともい

度維持など，農業への応用［１９］や，低温プラズマを利用した大気・水環境保全技術への応用［２０，２１］も期待されていう）の概念も共通である．真空装置内に存在する，例えば

る．

１０２０個の粒子（電子，イオン，中性の気体分子）のすべてを計算機内で取り扱うのは不可能であるので，例えば１０６個

２．
３ 低温プラズマシミュレーションの重要性の 超 粒 子 で 代 表 さ せ て 計 算 を 実 行 す る．こ の低温プラズマは非常に幅広い分野で実用化されている一１０２０/１０６＝１０１４，つまり一個の超粒子は１０１４個の実在粒子を方で，理解が及んでいない点も非常に多い．計測手段とし代表していることを表しており，これを超粒子の重みといては，比較的低圧のプラズマの場合，ある程度大きな空間

う．低温プラズマを取り扱う場合は，荷電粒子を取り扱うのプラズマが対象であれば，静電プローブや質量分析，

PIC/MCC 法の方が，中性粒子を扱う DSMC 法に比べて超LIF（レーザー誘起蛍光法）などが有効な場合もある．しか粒子の重みが小さくなるのが普通である．

しながら，近年応用を広げている大気圧プラズマの多く

１ステップの時間 !!は対象粒子種で大きく異なる．質は，その体積が非常に小さく，静電プローブのようにプラ量が軽く，加速によりエネルギーを得やすい電子の速度はズマに影響を与える測定方法は適さない．また．比較的実高いため PIC/MC 法の !!は１０−１０秒程度に短くとる必要が用的な測定手段の一つである発光分光測定についても，十ある．一方，DSMC 法の !!は１０−６秒程度である．

分な空間分解能を得ることは困難となる．更に，大気圧近移動と衝突を分離して取り扱う点は両者同様である．実傍の高圧力下では，プラズマによって様々な励起種やラジ際には，分子は !!の間のどこかで衝突をするのだが，計算カルが生成される上に，それらと他の原子・分子，イオン

では !!の時間，超粒子は衝突なしで移動し，!!後に（衝との反応によって二次的に生成される化学種の数が膨大と突すると判定された超粒子は）一斉に衝突するという取扱なる．例えば，N2 や O2 の振動励起状態を細かく分けて考えいをする．!!が平均衝突時間程度より小さければ，この取＋ないとしても，準安定な励起種が複数ある上に，O−
3 ，Ｏ4 ，

扱いは妥当である．

＋
N＋
3 ，N4 といったイオンも生成され，全ての化学種を定量

両者の衝突の取扱いは幾分異なる．PIC/MC 法ではバッ的に測定することは事実上不可能である．反応について

クグラウンドに存在する原料ガス（これは超粒子として取は，N++O2→NO++O のような異種のイオンと中性粒子のり扱わない）との衝突を考えるのに対し，DSMC 法では超反応の他，O+O2+M→O3+M（M は第三体粒子）のような粒 子 同 士 の 衝 突 を 取 り 扱 う．ま た，PIC/MC 法 で は 電

３体反応も無視できなくなる．粒子の平均自由行程は圧力

離，励起，解離などの反応を考慮する必要がある．最大のの上昇に反比例して短くなることから，反応はごく限られ

衝突周波数を見積もって仮衝突と実衝突の判定を行うな
２８８

Lecture Note

2. Fundamentals of Low-Temperature Plasma Simulation

K. Ikeda et al.

合である．ガスの圧力が (［Pa］で温度が 400 K のプラズマ

ど，計算効率をあげるための工夫には共通点がある．

!
"
(
!#
の場合，電子の平均自由行程は !
［m］と見積もられる

PIC/MC 法では電界を計算する必要があるという点もDSMC 法との大きな違いである．中性粒子の運動は電界のので［７］，系の代表長さを "［m］として !' &!
!
"という条件

影響を受けないが，荷電粒子の運動は電界の影響を考慮にに 代 入 す る と，("'!
!
#と な る．し た が っ て，例 え ば入れて運動方程式を計算することになる．

"$0.1 m の場合にはガス圧が 4 Pa 以上であれば電子に対一般的には粒子モデルに比べると流体モデルの計算負荷

して流体モデルを適用できそうである．逆に圧力がそれ以ははるかに小さいと言われているが［２６］，実はそうでもな下の場合は粒子モデルを用いたほうがよいであろう．ただい場合も多い．例えば，流体モデルには粒子種や反応式のし "の取り方には任意性があるので，これはあくまでも目数に比例して計算量が増える計算量の大きいルーチンがい安と考えてほしい．なおイオンの平均自由行程は電子よりくつかあるが，粒子モデルにはそれがない．粒子種や反応小さいので，同じ条件であれば流体モデルが適用できる．

式が非常に多い計算では，粒子モデルでの計算が速いケー

ところで，熱プラズマの場合には，電子とイオンと中性

スもある．また，流体モデルで計算する際の #)やメッシュ

粒子をまとめて一つの流体として取り扱えるので，運動方

幅も，粒子モデルに比べて常に大きくとれるというわけで

程式は一つだけ解けばよいが［３０］，非平衡低温プラズマのもない．

場合は粒子種ごとに異なった流体が粒子の衝突を介して相互に質量，運動量，およびエネルギーを交換しながら運動

PIC/MC 法と DSMC 法の詳細に関しては，紙面の都合する混合流体と考えなければならない．したがって，粒子上，いくつか文献を紹介するにとどめたい．
PIC/MC 法は１９９０年代初めには既に低温プラズマシミュ

種の数だけ方程式のセットを解かなければならない．しかレーションに用いられている［２８］．まずは Birdsall の教科

も現象が本質的に非定常であるため，定常解を一回求めて書［２２］を読むのがよいだろう．また，南部の著書［２７］に衝終わりという訳にはいかない．しかし，すべての方程式を突の取扱いについての記述がある．粒子モデルだけでなく忠実に解く場合には膨大な計算時間を要するため，実際に流体モデルにも当てはまるが，荷電粒子のシミュレーショ

は以下に示す近似を用いる場合が多い．

ンは非線形性が非常に強いため，空間メッシュ幅，時間ス

・電子とイオンについては，運動方程式を解くかわりにドテップ #)などの計算条件が適切に設定されていない場合

リフト−拡散モデルを用いて粒子フラックスを計算する．

には直後に計算が発散してしまうなどといった扱いづらい

・イオン温度の空間−時間変化は小さいとしてイオンに関するエネルギー方程式は解かない．
（イオンの温度は室面もある．

温などの既知であるとする．）

DSMC 法は１９６０年代前半に Bird により最初に開発され・中性粒子は全部まとめて一つの多成分気体とみなす．

た手法である．Bird は１
９７６年，１９９４年に DSMC 法の著書
［２４，２５］を出版している．また，２０１３年にも自費出版のサイトで新しい著書［２９］を発表している．南部の著書［２７］も取なお，荷電粒子の密度分布により電界が形成され，また電りかかりやすいであろう．DSMC 法は，原理的に計算が発界は外力の形で荷電粒子の運動に影響を与えるので，これ散することがないという流体シミュレーションにはない強らの方程式のほかに静電ポテンシャル（電位）に関するポ力な利点を持っている．ただし，妥当な計算結果を得るた

アソン方程式を同時に解く必要がある．電界は静電ポテン

めにはやはり，空間メッシュ幅や時間ステップ #)に注意

シャルの空間勾配から求められる．このような近似のもとが必要である．

で荷電粒子に関しては，電子および各イオン種についての連続の式と，電子に関するエネルギー方程式のみをポアソ

２．
５ 流体モデルによるシミュレーション手法

ン方程式とともに解けばよいことになる．

流体モデルでは，プラズマを連続体（流体）とみなし，質

ここで，ドリフト−拡散モデルについて説明しよう．例

量保存則を表す「連続の式」，運動量保存則を表す「運動方として電子を考える．外部磁界が無い（もしくは無視でき

程式」，そしてエネルギー保存則を表す「エネルギー方程る）場合，電子に関する運動方程式は以下のように書ける式」に基づいてシミュレーションを行う．未知変数は粒子（電子を表す下付の添え字は省略する）．

数密度 '$，流速ベクトル "$，温度 #$または圧力 ($であり，*"!""$!()&'!('!!'&&
'%$#*
!$
""#
(
)

これらは粒子種ごとに空間座標 +と時間 )の連続関数として定義される巨視的な量である．ここで添え字の $は粒子

（１）の 種 類 を 区 別 す る 指 標 で あ る．ま た 状 態 方 程 式

ここで，& は電子の質量，'は電子の数密度，"は電子の($ $'$%$#$，が成り立つと仮定する．%$ はボルツマン定数である．基礎方程式系は，これら未知変数の +および )速度，!は電界，$は素電荷，そして &は１個の電子が１秒に関する偏微分を含んだ非線形の偏微分方程式である．方ある．上式の左辺は電子流体の加速度を表し，右辺第１項

程式の具体的な形は例えば文献［３５］を参照されたい．このは圧力勾配による力，第２項は電気力，そして第３項は分ような連続体としての扱いが妥当であるのは，プラズマを子との衝突に起因する電子集団としての減速率をそれぞれ構成している粒子の平均自由行程 %とプラズマ容器サイズ

表している．イオンとの衝突による影響や粘性応力は無視

や電極間距離など系の代表的な長さ "の比 %""（これをク

している．いま，"が変化する時間的な尺度を 'としたとヌッセン数といい !%と記す）が１に比べて充分小さい場&%"であれば，左辺括弧内の第１項は右辺第３き，もし '

間に分子と弾性衝突をする平均回数（平均衝突周波数）で

２８９

Journal of Plasma and Fusion Research Vol.93, No.6 June 2017

項に比べて無視できる．また，)の大きさが熱運動の平均の粒子フラックスを求める際に，それを挟む両隣のセルで的な大きさ（!&&'程度）よりも小さければ，上式の右%""!

の密度の値を補間して用いる必要があるが，それには空間

辺第３項に比べて左辺括弧内第２項も併せて無視できる

１次元移流−拡散方程式の厳密解を利用する Scharfetter-［３１］．これらの場合には，式（１）から直ちに "が求められ

Gummel 法がよく用いられる［３４］．時間積分に関して陽解る．

法を用いる場合には，時間の経過とともに解が発散しない""! # !! ! %&
'%""'
&& &&
'

ように「電子が "(秒間にドリフト速度で進む距離がセルの幅を超えない」という条件が必要である．通常は，プロ（２）グラム内部で "(を自動的に決めるが，荷電粒子について

この式は，電子集団としての平均速度が局所的な電界に比は，おおよそ "("10−12〜10−10 秒程度である．中性粒子に例する速度（電界ドリフト速度）と，密度勾配に比例する

関しては，流れ場の変化する時間尺度が荷電粒子のそれと速度（拡散速度）の和で表されることを示している．式

比してオーダが大きく異なるため，荷電粒子の時間積分と（２）に密度を乗じてフラックス #（"'"）で表せば，は分離して，中性粒子計算のための別の "(を用いて計算#$ "#'$%$!!!$%'$を行う．これは，おおよそ "("10−8〜10−6 秒程度である．

（３）

荷電粒子を解く際には，中性粒子の密度分布を現在求めら&&を 移 動 度，と 表 す こ と が で き る．こ こ で，%$#!

れている最新の値に固定しておき，数周期分計算を進め!$%""!
&&を拡散係数という．ここでは，電子のみならる．その際，電子と分子の衝突によって生成される中性ラずイオンに対しても同様の近似を適用できるとして，粒子ジカル種などの生成率の空間分布（印加電圧１周期の平均種の区別を示す添え字 $を付けた．また，右辺第２項では値）を求めておく．続いて，中性粒子の計算を別の "(を用

温度を空間微分の外に出した形で近似している．右辺第１いて１［ms］程度の時間で計算を進める．その際には，今度項の複号は，正イオンに対してはプラス（＋），電子およびは荷電粒子種の密度の空間分布を固定しておき，また荷電負イオンに対してはマイナス（−）をとる．'に関しては，例

粒子の計算で求められた中性粒子の生成率の空間分布も併えば電極への 印 加電 圧 の 周 波 数 が $［Hz］で あ る ならばせて用いる．上記の通り，異なる時間スケールで行われる

1/ $［秒］と設定すればよい．式（３）において，!は瞬時の荷電粒子の計算および中性粒子の計算のやり取りを反復す&! 1 の場合においても，緩やかに変化す電界であるが，'ることによって（周期的）定常状態時の低温プラズマの諸る実効電界 !#$$を !の代わりに用いることで，式（３）と同特性を求める．

じ形式で粒子フラックスを計算する手法が提案されている［３１］．このとき，移動度や拡散係数および各種の衝突周波

参考文献

数などを輸送パラメータや輸送係数，スオームパラメータ［１］K. Ikeda et al., J. Vac. Soc. Jpn. 50, 423 (2007).
［２］E. Pfender, Gaseous Electronic, Vol. 1, M.N. Hirsh and H.
J. Oskam, Eds. (New York: Academic Press, 1978) p.302.
［３］F. Tochikubo et al., J. Plasma Fusion Res. 92, 680 (2016).
［４］Y. Tanaka et al., J. Plasma Fusion Res. 82, 479 (2006).
［５］Y. Uesugi, J. Plasma Fusion Res. 89, 207 (2013).
［６］F. Tochikubo et al., J. Plasma Fusion Res. 91, 307 (2015).
［７］菅 井 秀 郎：プ ラ ズ マ エ レ ク ト ロ ニ ク ス（オ ー ム 社，
２
０
０
０）
．
［８］大竹尚登：DLC の応用技術 ‐進化するダイヤモンドラ
イクカーボンの産業応用と未来技術‐（シーエムシー，２０１
６）．
［９］http://www.spstj.jp/publication/archive/vol19/Vol19
sub size 6 down 10 No6 sub size 6 down 10 1.pdf.
［１０］松崎紘一 他：日口腔インプラント誌 27, 30 (2014).
［１１］是津信行 他：表面技術 65, 227 (2014).
［１２］津山青史：鉄と鋼 100, 71 (2014).
［１３］T. Segawa et al., J. Plasma Fusion Res. 91, 665 (2015).
［１４］H. Nishida et al., J. Plasma Fusion Res. 91, 661 (2015).
［１５］http://plasmamedizin.com/［１６］S. Ikehara et al., Plasma Process Polym. 12, 1348 (2015).
［１７］K. Miyamoto et al., J. Clin. Biochem. Nutr. 60, 25 (2017).
［１８］T. Hirata et al., J. Plasma Fusion Res. 91, 771 (2015) .
［１９］K. Takagi et al., J. Plasma Fusion Res. 90, 531 (2014).
［２０］M. Okubo et al., J. Plasma Fusion Res. 84, 121 (2008)［２
１］A. Yamatake et al., IEEE Trans. Plasma Sci. 34, 1375などという．これらは，原料ガスの種類と換算電界（電界強度をガス分子の数密度で除した値）に応じて異なるため，式（３）において一般に定数ではない．電子はイオンに比べて軽いため，スオームパラメータを理論的に求める方法は両者で異なる．これに関しての詳細は，電子に関しては文献［３２］を，またイオンに関しては文献［３３］をそれぞれ参照されたい．いずれにしても，スオームパラメータを求めるためには，低温プラズマシミュレーションで考慮する荷電粒子と中性分子（原子）との各種衝突断面積（一般に衝突する２粒子の相対速度に依存する）のデータをあらかじめ準備しておく必要がある．またスオームパラメータが実験的に求められている場合には，その値を利用することもある．
実際に低温プラズマシミュレーションを行うには，微分方程式を離散化しなければならない．空間に関してはメッ
シュ分割を行って，計算対象となる領域を多数の（小さい）空間領域（セル）に分割する．未知変数は各セルにおいて一つの値を持つことになる．差分法や有限体積法を用いて方程式の離散化を行う場合，密度や温度などのスカラー量
をセルの中心で定義し，粒子フラックスなどのベクトル量
は，隣あうセルとの境界面においてその面に垂直な成分で定義する場合が多い［３４］．このように諸量を設定すると，(2006).
［２２］C.K. Birdsall and A.B. Langdon, Plasma Physics via Com-

セル内のスカラー変数の時間変化が，セル表面を通過する
粒子フラックスで表現できるため都合が良い．セル表面で
２９０

Lecture Note

2. Fundamentals of Low-Temperature Plasma Simulation

puter Simulation (McGraw-Hill, 1985).
［２３］K. Nanbu, IEEE Trans. Plasma Sci. 28, 971 (2000).
［２４］G.A. Bird, Molecular Gas Dynamics (Clarendon Press,

K. Ikeda et al.

いけ

けい

圭

１
９
８
９年東京理科大学卒業，同年日電アネル
バ㈱（現，キヤノンアネルバ㈱）入社．２
００１
年ウェーブフロント入社．２
０
０
８年１
０月㈱ア
テナシス設立．化学反応を含んだマルチ
フィジックス解析に従事．今回，小田先生より貴重な機会をいただきました．これから低温プラズマシミュレーションを始めようとしている方々に少しでも参考にしていただければ幸いです．

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

だ池 田!!!!!!!!!!!!!!!!!!!!!!!!!!

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

まつ

なが

ふみ

ひこ

１
９
８
２年京都大学理学部卒（おもに数学を専攻）．同年セン チ ュ リ リ サ ー チ セ ン タ ㈱（現伊藤忠テクノソリューションズ㈱）に入社し，おもに流体に関するシミュレー
ションを担当する．２
０
０
２年にペガサスソフトウェア㈱に移り，流体モデルを用いた非平衡プラズマシミュレーション用のソフトウェア開発に従事し現在に至る．

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

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

松 永 史 彦

だて

ひろし

中 舘

なか

博

１
９
９
８年 東北大学大学院工学研究科 原子核工学専攻 後期課程 修了．博士（工学）．株
式会社 CRC を経てペガサスソフトウェア
株式会社を設立．４
０代後半にして社内最年
少．粒子モデルのシミュレーションが専門．GUI やインス
トーラ回りだけでなく，会社のウェブ更新，ゴミ出し，電気
ポットにお湯の補給などの雑用も専門．

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

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

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

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

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

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

1976).
［２５］G.A. Bird, Molecular Gas Dynamics and the Direct Simulationof Gas Flow (Clarendon Press, 1994).
［２６］南部健一：プラズマ・核融合学会誌 77, 1137 (2001).
［２７］南部健一：原子・分子モデルを用いる数値シミュレー
シ ョ ン 第３章 モ ン テ カ ル ロ 法 の 基 礎（コ ロ ナ 社，
１
９
９
６）
．
［２
８］C.K. Birdsall, IEEE Trans. Plasma Sci. 19, 65 (1991)［２９］G.A. Bird, The DSMC Method, createspac, https://www.
createspace.com/3689652 (2013).
［３０］Watanabe et al., Plasma Fusion Res. 87, 514 (2011).
［３１］Gogolides et al., J. Appl. Phys. 72, 3971 (1992).
［３２］G.J.M. Hagelarr et al., Plasma Sources Sci. Technol. 14,722 (2005).
［３３］Mason et al., Transport Properties of Ions in Gases (Wiley,1992).
［３４］真壁利明：プラズマエレクトロニクス（培風館，１
９
９
９）
．
［３
５］Lieberman et al., Principles of Plasma Discharges and Materials Processing (Second edition, Wiley, 2005).

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

２９１

J. Plasma Fusion Res. Vol.93, No.7 (2017)３
２３‐３２６

講座 実践低温プラズマシミュレーション

３．低温プラズマのシミュレーション事例
3. Simulation Examples of Low-Temperature Plasma

３．
１ 粒子モデルによるプロセスプラズマシミュレーション
3.1 Processing Plasma Simulation Using Particle Model
田中正明
TANAKA Masaaki
ペガサスソフトウェア株式会社（原稿受付：２
０
１
７年５月１
１日）

粒子モデルを使用したプロセスプラズマシミュレーションの適用事例として，プラズマを利用する物理蒸着の一つであるマグネトロンスパッタリングにおけるプラズマ，表面改質のためのプラズマイオン注入成膜法におけるプラズマ，摩擦現象における滑り接触で生成されるトライボプラズマを紹介する．粒子モデルによるプロセ
スプラズマシミュレーションは，装置サイズからサブミクロンサイズまで，また衝突プラズマから無衝突プラズ
マまで，その使用範囲は広く，多くの産業分野で利用されている．
Keywords:
plasma simulation, particle model, magnetron sputtering, diamond-like carbon, triboplasma

３．
１．
１ はじめにmethod）［２‐５］が，また外力項が無い中性粒子についてはプラズマの主な応用分野として，核融合並びに多くの産

DSMC 法（Direct Simulation Monte Carlo method）［５，６］が業分野で利用されているプラズマプロセスがある．プラズ

使用される．右辺の衝突項については，中性粒子間衝突，マプロセスとは，プラズマを用いた材料の加工（成膜，

荷電粒子−中性粒子間衝突，クーロン衝突について詳細なエッチング，表面改質）のことである．材料の加工におい記述がある［５］．

ては，ガスの温度は常温に近く，プラズマとしては電離度

上記粒子シミュレーションは，ボルツマン方程式の未知が低く（数％程度），電子温度は高いが，中性粒子の温度が量である速度分布関数を求め，種々の物理量は得られた速低いことが求められる．このようなプラズマは非平衡・低度分布関数から求められる．したがって連続体モデルにお温プラズマと呼ばれ，平衡・高温プラズマである核融合プ

ける速度分布の仮定，および輸送パラメータが不要であ

ラズマとは性質が大きく異なる．非平衡・低温プラズマ

る．また粒子シミュレーションは，境界条件の取り扱いがは，材料加工の他，半導体製造，医療，バイオ，環境，人容易である．しかし，プラズマシミュレーションに不可欠

工衛星推進などの分野で拡大している．本節ではこの非平な気相反応，表面反応に関するデータが重要であるが，薄衡・低温プラズマをプロセスプラズマという．低温プラズ

膜生成のためのオルトケイ酸テトラエチル（Si（OC2H5）4；

マについては，本学会誌の小特集「プラズマ応用技術におTEOS），テトラメチルシラン（Si（CH3）4；TMS），ヘキサ

けるシミュレーション研究」の中の「２．低温プラズマ放

メチルジシロキサン（（TMS）
2O；HMDSO）などの高分子

電の基礎」［１］を参照されたい．流れの代表長さと分子の平についてはシミュレーション可能な気相反応データが，ま均自由行程の比で表されるクヌッセン数が 0.01 以上となる

た表面反応についても表面の状態により反応が大きく変化

中間流，分子流においては粒子間衝突が減少し，非平衡性するため，統一的に適用可能なデータは存在しない．

が強くなり，連続体とみなされずボルツマン方程式に立ち

３．
１．
２ マグネトロンスパッタ装置への利用例

返って中性粒子・荷電粒子の運動を扱うことになる．このようなプラズマの支配方程式は，ボルツマン方程式とマクマグネトロンスパッタ法では，ターゲット表面に磁界をスウェル方程式である．そしてその数値解法として代表的印加して，イオン入射による二次電子をローレンツ力で捉な粒子シミュレーション手法は，ボルツマン方程式の左辺

えてサイクロイドまたはトロコイド運動させることによりの項である外力項として電磁界を考慮する荷電粒子についアルゴンガスとのイオン化衝突の頻度を増大させ，ターては，PIC/MC 法（Particle-In-Cell/Monte Carlo collision

ゲット付近に高密度プラズマを生成させることで，成膜速PEGASUS Software Inc.
author’s e-mail: mtanaka@psinc.co.jp

３２３!2017 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.93, No.7 July 2017

度の高速化を可能にする．また，電子が磁界による束縛から逃れて基板に入射するまでにイオン化衝突による運動エ
ネルギー低減が十分なされるために，高エネルギー電子の基板衝撃は起こらずそれに伴う基板温度上昇が抑制される．これによって，マグネトロンスパッタ装置は急速に普及したが，膜厚・膜質の均一化などに課題があることから，さまざまなマグネトロンの方式が開発されている．このうち，現在ではスパッタリング装置においても枚葉式の装置が主流となっており，プレーナ型が最も普及している．主に金属ターゲットにおける直流マグネトロン放電は平板上陰極面に平行に磁場をかけて放電させる．低い圧力

（a）
解析モデル

（1 Pa 以下）でも，電子の "!!ドリフトによる周回運動の効果（マグネトロン効果）によりドーナツ状の高密度プラズマが生成される．電流密度が高く，イオンは高エネル
ギー（数百 eV）で陰極をたたき，スパッタリング粒子の平均自由行程が長いため薄膜形成などに使われている．また
反応性スパッタリング，化合物ターゲットの場合は直流放電ではなく，RF 放電が利用される．近年スパッタリングの一種で，低い duty 比（１パルス周期の中でパルスが"オン"になっている時間の割合）で集中させた高い電力を，瞬間（b）電子密度分布（XY 断面）

的にカソードに投入することで高密度のプラズマを形成することにより 表 面の 平 坦 性 が 良 く，緻 密 な 膜 が できる
HiPIMS（High Power Impulse Magnetron Sputtering）［７］の研究が盛んに行われている．マグネトロンスパッタシ
ミュレーションでは，ターゲットへのイオン入射による二次電子放出係数，スパッタ率，スパッタ粒子の放出エネル
ギー分布，放出角度分布を与えることは特に化合物のス
パッタリングにおいて難しく，イオン−固体相互作用の研究は今後も続くであろう．ここではアルゴンガス圧 0.4 Pa，
ターゲット印加電圧−300 V，最大磁場強度 350 G の直流マ
グネトロンプラズマの３次元シミュレーション例として，（c）
電子密度分布（XZ 断面）

（d）
電位分布

解析モデル，電子密度の XY 断面分布，XZ 断面分布，そし図１て電位の XY 断面分布を図１（a）〜（d）にそれぞれ示す．

直流マグネトロンプラズマの空間３次元シミュレーション
事例．

３．
１．
３ プラズマイオン注入成膜法への利用例
DLC（Diamond-Like Carbon）膜は，表面保護としての耐Ion Implantation and Deposition; PBII&D 法）のシミュレー

摩耗性，低摩擦係数，耐腐食性（機械，自動車部品，加工

ションを紹介する．PBII&D には，RF パルスと DC パルス

用ドリル，HDD ディスク＆ヘッド），酸素透過防止としてを組み合わせた方式［１０］と，正・負の DC パルスを使用すのガスバリア性（ペットボトル），有害物質の抑制としてるバイポーラパルス方式［１１‐１３］がある．そして各パルスの生体親和性（ステント），赤外線透過性（赤外線窓）など

間のアフターグローの時間は数十 μs である．このシミュに利用されている．特長として，（１）室温近傍の低温で成レーションでは，バイポーラパルス方式を採用し，正パル膜が可能，（２）基板の種類が多く非耐熱性のプラスティッスとして 1 kV を，そして負パルスとして−3 kV をそれぞ

クにも成膜が可能，（３）大面積かつ複雑な３次元形状にもれ 1.5 μs 印加した．その際，荷電粒子２７種類，中性粒子２４種成膜が可能，（４）成膜しても表面粗さが変化しない，など

類，気相反応７０９種類をそれぞれ考慮した．また，初期密度が挙げられる．

14−3は電子および C2H＋とした．正パルス終
2 共に 1.0×10 m

DLC 成膜法として，イオンビーム蒸着法，アークイオン

了 後 の 電 子 密 度，C2H＋
2 密 度 の 最 大 値 は そ れ ぞ れ

プレーティング法，マグネトロンスパッタリング（アンバ

9.0×1014 m−3，7.9×1014 m−3，電極への入射イオンフラッ

ランス型）法，プラズマ CVD 法，そしてプラズマイオン注

クスは入射面平均で 1.7×1018 m−2s−1 である．負パルス終入成膜法などがある．また近年マイクロ波プラズマによる了 後 の 電 子 密 度，C2H＋
2 密 度 の 最 大 値 は そ れ ぞ れ

高速成膜法も研究されており［８］，DLC 成膜プロセスにつ

7.8×1014 m−3，6.9×1014 m−3，電極への入射イオンフラッ

いて詳細な記述がある［９］．

クスは入射面平均で 1.6×1019 m−2s−1 である．図２（a）〜（g）に，解析モデル，正パルス終了後の電 子 密 度，C2H＋2 密

本節では，プラズマイオン注入成膜法（Plasma Based
３２４

Lecture Note

3.1 Processing Plasma Simulation Using Particle Model

M. Tanaka

３．
１．
４ トライボプラズマへの利用例
１９５０年代の後半に，英国のケンブリッジ大学で摩擦に伴う閃光温度が計測されて以来，トライボロジーの諸問題は摩擦に伴う温度上昇をもとに解析され，解決されてきた．
しかしながら，熱エネルギーのみによっては説明できない現象が多々観察され，トライボロジーの技術開発を妨げてきた．その代表的なものが説明のできないトライボケミカ
ル反応である．このことは摩擦の熱エネルギー以外に高いエネルギー状態が存在すことを示唆していた．一方，固体の新生面からエキソ電子と呼ばれる電子が放出されること（a）解析モデル（b）電子密度分布（正パルス終了後）が知られており，このことは摩耗新生面から電子が放出されることを示唆していた．そこでこの電子に着目し，大気中，油中で計測可能な摩擦に伴って放出される電子，次いで負の荷電粒子，正の荷電粒子，さらに光子を計測することのできる新奇なトライボエミッション計測装置を用いて，これらの素粒子や荷電粒子が摩擦に伴って放出されることが明らかになった．この摩擦にともなう粒子放出現象
をトライボエミッションと呼ぶ．正および負の荷電粒子を同数検出したことから，摩擦接触点にプラズマ（トライボ
プラズマ）が発生していると結論づけられた．このトライ＋（c）C2H2 密度分布
（正パルス終了後）

ボプラズマは摩擦帯電による高電界により周囲気体の放電（d）電位分布（正パルス終了後）によって摩擦接触点の隙間に発生することを突き止め，その撮像に成功し，これを発見している［１４］．このトライボ
プラズマこそトライボケミカル反応を始めとする様々なト
ライボロジー問題を引き起こす主要な原因の一つと考えられるに至った．トライボプラズマのトライボロジー問題解決に向けた応用技術の展開が期待されている．しかしながら，トライボプラズマの応用技術を開発するには，トライ
ボプラズマを構成するこれらの諸粒子の発生強度や発生分
布を明らかにすることが不可欠であるが，実験的には極め（e）電子密度分布（負パルス終了後）て困難である．一方，PIC/MC シミュレーション法はこれ＋（f）C2H2 密度分布
（負パルス終了後）

らの諸粒子の発生強度や発生分布を明らかにすることのできる極めて有力な武器であり，本法を適用してトライボプ
ラズマの解析を試みた．トライボプラズマの発生分布が一致しPIC/MC法がトライボプラズマの解析や応用技術の展開に有用であることを示すことができた［１５］．その後，
データベースを再構築し，荷電粒子１９種類，中性粒子１９種類，気相反応３５８種類で再計算を行い，発生強度分布が実測と一致することを確認した．荷電粒子密度分布の一部とし＋＋＋て，図３（a）〜（g）に解析モデル，電子，N＋
2 ，N4 ，O2 ，O ，O−，の密度分布をそれぞれ示す．これら荷電粒子密度の 最 大 値 は そ れ ぞ れ 8.9×1021 m−3，1.1×1021 m−3，（g）電位分布（負パルス終了後）
図２

4.2×1021 m−3 ，1.4×1021 m−3 ，1.4×1021 m−3 ，6.9 × 1019 m−3であった．これらの図から電子密度が最大となるギャップ

プラズマイオン注入成膜法（PBII&D 法）のシミュレーショ
ン事例．

長はパッシェン則から得られる値とほぼ同値であり，バルクプラズマは上記の正負の荷電粒子でほぼバランスしてい度，電位の各空間分布，そして負パルス終了後の電子密度，

ることがわかる．

C2H＋2 密度，電位の各空間分布をそれぞれ示す．これらの図

一方，油潤滑下においても，トライボプラズマの撮像にから，正パルスにより被覆物近傍に高密度プラズマが生成成功し，油潤滑下でもトライボプラズマが発生することをされ，その後の負パルスにより被覆物にイオンが均一に入証明し，その発生メカニズムを明らかにしている［１６］．実射していることがわかる．

用上とりわけ重要なのは，この油潤滑下で発生するトライ
ボプラズマであり，このトライボプラズマの油剤への作用
３２５

Journal of Plasma and Fusion Research Vol.93, No.7 July 2017

トライボプラズマの観点からさまざまなトライボロジーの諸問題が解決されるものと期待される．

３．
１．
５ まとめと今後の課題
近年多くのコアを搭載した比較的安価な PC，ワークス
テーションが販売されていることから，並列化効率の良い粒子モデルがよく使用される．粒子モデルによる計算時間
（a）
解析モデル

は，扱う粒子種数にはそれほど依存せず，超粒子（サンプ（b）
電子密度分布ル粒子）の総和に依存する．荷電粒子の粒子モデルではマクスウェル方程式もしくはその一部であるポアソン方程式を解く計算時間をいかに少なくするかによる．またプロセ
スプラズマでは様々な分野で利用されるため，電子，イオ＋

ンを流体モデルと粒子モデルに分けるハイブリッドモデル＋（c）
N2 密度分布（d）
N4 密度分布が今後さらに必要となる．そして原子・分子データベースの構築・整備も必要不可欠である．

参考文献＋［１］伊達広行：プラズマ・核融合学会誌 80, 113 (2004).
［２］C.K. Birdsall and A.B. Langdon, Plasma Physics via Computer Simulation (McGraw-Hill, 1985).
［３］C.K. Birdsall, IEEE Trans. Plasma Sci. 19, 65 (1991).
［４］R.W. Hockney and J.W. Eastwood, Computer Simulation
Using Particles (Taylor and Francis, Inc., 1988).
［５］K. Nanbu, IEEE Trans. Plasma Sci. 28, 971 (2000).
［６］G.A. Bird, Molecular Gas Dynamics and the Direct Simulationof Gas Flow (Clarendon Press, 1994).
［７］小島啓安：現場のスパッタリング薄膜 Q&A 第２版（日刊工業新聞社，２０
１５）．
［８］H. Kousaka et al., IEEE Trans. Plasma Sci. 41, 1830 (2013).
［９］上坂裕之：プラズマ・核融合学会誌 90, 76 (2014).
［１０］http://pekuris.co.jp/jyoho̲pbii.html
特開２
０
０
１
‐
０
２
６
８
８
７，２
０
０
１
‐
２
０
７
２
５
９，２
０
０
４
‐
３
２
３
９
７
３
［１
１］宮川佳子：プラズマ・核融合学会誌 80, 126 (2004).
［１２］W. Park et al., J. Phys. D : Appl. Phys. 47, 335306 (2014).
［１３］Y. Hirata et al., Int. J. Refract. Met. Hard Mat. 49, 392 (2015).
［１４］K. Nakayama and R.A. Nevshupa, J. Phys. D: Appl. Phys.
35, L53-6 (2002).
［１５］K. Nakayama and M. Tanaka, J. Phys. D: Appl. Phys. 45,
495203 (2012).
［１６］K. Nakayama, Tribol. Lett. 41, 345 (2011).

＋（e）
O2 密度分布（f）
O 密度分布（g）
O−密度分布
図３

摩擦現象における滑り接触で生成されるトライボプラズマのシミュレーション事例．

により引き起こされるトライボケミカル反応であって，電子やイオン，ラジカルなどのトライボプラズマ中の活性中間体を解析することが不可欠である．
そ こ で，潤 滑 油 の 基 油 を 構 成 す る 飽 和 炭 化 水 素 分子
（CnH2n＋2）を CH4 ガスでシミュレーションし，電子，様々なラジカル，イオン，反応生成物の発生強度や発生分を明らかにすることができた．このシミュレーションで水素の発生量が最も大きく現れたので，これに基づき実際に炭化水素油中で放電させたところ水素が検出された．本結果から，このシミュレーションによるトライプラズマ粒子解析
を通じてトライボロジー問題解決に有用であることが示された．PIC/MC シミュレーションと合いまった実験により!!!!!!!!!!!!!!!!

た

なか

まさ

あき

田 中 正 明!!!!!!!!!!!!!!!!

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

１
９
７
４年東京理科大学理学部卒．同年セン
チュリリサーチセンタ㈱（現伊藤忠テクノ
ソリューションズ㈱）に入社，数値解析，
統計解析ソフトウェア開発，核融合炉ダイ

バータプラズマシミュレーションに従事する．２
０
０
２年にペガ
サスソフトウェア㈱を設立し，希薄気体，非平衡プラズマ，
形状シミュレーション等のソフトウェア開発に従事し現在に至る．

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

３２６

J. Plasma Fusion Res. Vol.93, No.7 (2017)３
２７‐３３１講座 実践低温プラズマシミュレーション

３．低温プラズマのシミュレーション事例
3. Simulation Examples of Low-Temperature Plasma

３．
２ 流体モデルによるプロセスプラズマシミュレーション
3.2 Processing Plasma Simulation Using Fluid Model池田圭

IKEDA Kei
株式会社アテナシス（原稿受付：２
０
１
７年５月１
１日）

流体モデルを利用したプラズマシミュレーションの適用事例として，半導体製造プロセスへの応用例を２例紹介する．一つは，半導体製造装置でしばしば利用されるガス分散板（シャワーヘッド）の孔に着目し，その孔の出口近傍でプラズマが局所的に増幅されていることを見出し，シャワーヘッドが損傷する原因の究明へと繋がった例である．もう一つは，プラズマ CVD プロセスへの応用例であり，従来知られているプロセスガスに Arを加えると同時に，ステージに接続されたインピーダンスを調整することによって，基板に入射するエネルギーを制御し，堆積する膜のステップカバレッジの向上と低抵抗化の両方を同時に達成した．どちらの事例もプラズ
マシミュレーションが製品開発への指針となった実用的な例である．
Keywords:
plasma simulation, fluid model, semiconductor, etching, chemical vapor deposition, showerhead

３．
２．
１ はじめにロセスは流体モデルで扱うことが可能であり，計算コスト

半導体製造プロセスは，デバイスを作成する過程で必要も比較的安価である．また，ハードウェアの性能向上と並となる微細加工の高精度化と基板サイズの大型化にともな

列化の利用により，近年は取り扱える化学種や反応式の数

い，従来構造の延長上にある装置開発から，シミュレーも増え，装置開発の指針となり得るレベルになってきた．

ションを活用したより効率的な装置開発へとシフトしてき

本節では，半導体製造装置開発に流体モデルによるプロている．解析対象も多岐に渡り，熱流体解析は勿論，電磁セスプラズマシミュレーションを利用した Denpoh［４，８］

場や反応を考慮したプロセスシミュレーション等，流体モによる２例について主に紹介する．一つはエッチング装置

デルを利用した様々な検討が行われている［１‐８］．ドライに応用した例［４］で，シャワーヘッドのガス孔近傍に着目

エッチングと呼ばれる微細加工技術は，基板表面からエッし，プラズマがガス孔の出口近傍で局所的に増幅されるこ

チングされた副生成物がプラズマ中に入り，そのまま戻るとを明らかにしている．従来の構造では，シャワーヘッド

ものもあれば，反応して別の化学種となって再び基板に付の孔の出口付近が損傷し，シャワーヘッドの寿命を短くす

着する化学種もあることから，その気相および表面反応モる問題があったが，その原因が解明された．もう一つは，

デルは極めて複雑であり，その構築も容易ではない．エッ成膜プロセスの一つとして重要なプラズマ CVD（Plasma-

チングレートは基板に入射するイオンエネルギーや温度に

enhanced Chemical Vapor Deposition）装置に応用した例

対する依存性を有するが，汎用的な表面反応モデルの構築［８］で，シミュレーションを活用することにより，Denpohに必要な実験データは現在もほとんど揃っていない．そのの共著らと共に，ステージに接続されたインピーダンスをため，プラズマからの入射粒子と基板表面の原子・分子の調整することによって基板に入射するイオンの密度とエネ間で生じる反応を分子動力学（Molecular Dynamics : MD）

ルギーを制御する方法を確立し，堆積する Ti のステップカを用いて検討する試みも進められている［９］．しかしながバレッジと膜の低抵抗化の両方を同時に達成することに成ら，MD シミュレーションは計算コストが非常に高く，現功している．

在も基礎的な研究が進められている段階である．装置シ

３．
２．
２ エッチング装置への応用例

ミュレーションに関しては，前章で述べた通り，大きく分けて粒子モデルと流体モデルの２通りがあるが，多くのプエッチング装置の多くは容量性結合型プラズマ（Capaci-

ATHENASYS Co.,Ltd.
authors’ e-mail: dk@athenasys.co.jp

３２７!2017 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.93, No.7 July 2017

tively Coupled Plasma ; CCP）と誘導性結合プラズマ（In-なる［４］．圧力は一般的なプロセスよりも低い 25 mTorr

ductivity Coupled Plasma ; ICP）に区分される．一般に ICP（3.33 Pa）
，シ ャ ワ ー ヘ ッ ド に 印 加 さ れ る 電 圧 はは数Pa程度の低圧力に用いられることが多い一方で，CCP##$"$"#!"(%)&'"!
#に 設 定 さ れ て い る．こ こ で，はおよそ 10〜1000 Pa のやや高い圧力で用いられる［１０］．

"(%!150 V，周 波 数 #は 60 MHz で あ る．こ の と き，Ar

プラズマの生成方法が大きく異なり，特にフラットパネル

ガスが孔の上部から供給され，中心軸と対向する外側（電のような大面積の基板処理には，CCP タイプの装置が多く極間）の境界から排気される．

用いられる．CCP 装置の基本的な構造は平行平板式となっ本シミュレーショ ン に は 有 限 体 積 法（Finite Volumeており，１つもしくは複数の基板がステージの上に配置さMethod ; FVM）ベースとするマルチフィジックスソフト

れ，上部には基板に対抗するように多数の小さな孔を有すウェア CFD-ACE＋［１２］が用いられており，流体モデルをるシャワーヘッドが配置される．導入するプロセスガスはベースとした解析が可能である（新しいバージョンでは，

このシャワーヘッドを介して供給され，同時にシャワー

イオンのエネルギーおよび角度分布について，粒子法を用ヘッドが電極の一つを兼ねている．電圧を印加する従来のいてポスト解析する機能も搭載されているが，本計算が実方法は，シャワーヘッド側に印加されるアノードカップル

行された時点ではその機能は搭載されていない）．プラズ

型とステージ側に印加されるカソードカップル型の２通り

マと中性粒子の流れはカップルされており，電位の計算にであったが，近年は両電極に印加されるタイプが増えていはポアソン方程式を解いて求めている．なお，電子の輸送

る［６］．シャワーヘッドに 13.56，27 も し く は 60 MHz と係数と反応レートについては SWARM が用いられている

いった工業用の高周波（Radio Frequency: RF）を印加する［４］．一般的に，プラズマシミュレーションを報告した例

ことによりプラズマを生成・維持し，ステージ側には主にでは，ガス密度を一定として計算しているものも少なくな

イオンエネルギーを制御する目的で，イオンが追随し始め

いが，本シミュレーションは計算領域内の圧力や密度の変る 1 MHz 程度より低い周波数［１１］を印加する２周波印加化が重要であるため，別途状態方程式を考慮している．な

CCP も検討されている．片方の電極に異なる周波数を重畳

お，CCP タイプの装置では，一般にプラズマの体積 " とプする方式もあり，いずれも入射する荷電粒子のタイミング

ラズマが接する壁の面積 !の比"!
!が小さく，壁からの熱やそのエネルギー・ラジカルの密度等を制御する方法としのロスが大きいために，ガス温度の上昇は小さいことが期て用いられている．Denoph の報告例は，基板側を接地し，待される．そのため，本シミュレーションではガス温度一上部のシャワーヘッドに RF が印加されるアノードカップ

定として計算されている．以下に Denpoh による文献［４］でル型の構造となっている．

示されている代表的な計算結果を示す．

計算モデルはシャワーヘッドの孔の一つに着目したもの

孔を流れるガス流量 0.9 sccm の条件において，1/4 周期

で，基板とシャワーヘッドの間隔は 35 mm，"0.5 mm の吹の位相におけるイオン化レートを図１に示す（ガス流量

き出し孔を中心軸に配置した二次元軸対称モデルである．

0 sccm の条件および 0.9 sccm の他の位相におけ る 結 果

以下で説明する計算結果を示した図１の片側が計算領域とは，Denpoh の文献［４］を参照）．計算モデルは軸対称であるが，反対側にミラーした結果が示されている．ガス流量
0.9 sccm の時には，孔の出口で局所的にイオン化レートが大きくなったが，ガスが流れていない条件では，全ての位相で局所的なイオン化レートの増幅は見られていない．
次に，中心軸に沿った圧力分布を図２に示す．点線の位

図１

図２

イオン化レート（1/4 周期，ガ ス 流 量 0.9 sccm）．（Reprinted from “Locally enhanced discharges at gas holeoutlets of a showerhead in a plasma etching reactor”,
K. Denpoh 2009 J. Phys. D: Appl. Phys 42, 032003 (2009),! IOP Publishing. Reproduced with permission. All rights
reserved）

３２８

中 心 軸 上 の 圧 力 分 布（孔 の 出 口 を 横 軸０と し て 表 示）．
（Reprinted from “Locally enhanced discharges at gas
hole outlets of a showerhead in a plasma etching reactor”,
K. Denpoh 2009 J. Phys. D: Appl. Phys 42, 032003 (2009),! IOP Publishing. Reproduced with permission. All rights
reserved）

Lecture Note

3.2 Processing Plasma Simulation Using Fluid Model

K. Ikeda

置よりも左側が孔の内部（上流側）で，点線よりも右側がよびジュール加熱を加えた計算が必要となる．更に，表面電極間を示す．孔の内部では，コンダクタンスが非常に小反応におけるイオンや活性種のエネルギー損失も考慮するさいために大きな圧力勾配が生じている．本シミュレー

ことにより，より高精度の解析モデルとすることも可能で

シ ョ ン の 結 果，孔 の 出 口 付 近 に お い て も 約 700 mTorr

ある．Miyashita ら［７］は，CFD-ACE＋を用いて石英製の（93.3 Pa）まで圧力が上昇しており，プラズマのバルクの圧真空容器を用いた ICP 装置の解析を行い，コイルに掛かる力と比較してはるかに高いため，イオン化レートも増幅さ電位の影響を考慮して CCP のタイプと同様にプラズマ電

れて電子密度もその近傍で局所的に高くなっていることが位を求め，シースで加速されるイオンのエネルギー分布・

示されている．

角度分布を利用して石英がスパッタエッチされる分布を計

中心軸上に沿った電子とイオンの数密度分布を図３に示算した．実験結果から得られた石英の損傷の様子（スパッす（横軸は図２と同様である）．位相は 1/4 周期の時点であタエッチレートの分布）とシミュレーション結果を比較しり，実線は電子密度，点線はイオン密度を示す．ガスが導たところ，その傾向が良く一致し，石英が損傷する分布の入されている場合，バルクのプラズマ密度も全体的に上昇メカニズムを明らかにしている．このように，装置の部材

していることがわかる．これは，位相が進むにつれてシャが損傷する原因を探り，その対策を検討することは装置設

ワーヘッドがアノードからカソードへと相対的に変わる際

計の上で重要な課題であり，現在のシミュレーションは有に，孔の出口近傍で生成した二次電子がバルクに戻り，プ用なツールとなっている．

ラズマ密度の上昇に寄与したことを示唆する．

３．
２．
３ プラズマ CVD 装置への応用例

文献［４］では，実際の装置で実験を行い，孔の断面形状を比較した結果も示されている．ガスを流した条件（a）では

Ti 薄膜はバリア膜やシード層，密着性向上のための接着ガスを流さない条件（b）と比べ，孔の出口近傍が広がって層等に利用されており，その堆積方法の一つとしてプラズ

おり，図３の結果と対応することが確認されている．このマ CVD が利用されている．熱 CVD を利用した製膜も可能

現象は，ステージ側に高周波を印加したカソードカップル

であるが，プラズマを利用することで低温化が可能である型でも確認され，電圧の印加方法には依存しなかったこと

ことも低温プラズマを利用する利点の一つである［１４］．プも報告されている．

ラズマ CVD による Ti 製膜プロセスは，TiCl4/H2 をベース

以上のように，エッチング装置におけるシャワーヘッドとしたガス系が以前より知られている［８］．これは，気相の孔の出口近傍では局所的にプラズマの密度が上昇し，バ中で H2 から解 離・生 成した H が TiCl4＋H→TiCl3 + HCl，

ルクのプラズマ密度の増加にも寄与していることがシミュ

TiCl3 + H→TiCl2 + HCl のように反応し，更に TiClx が表面

レーションによって明らかにされた．エッチング装置で利に吸着して TiCl（s）＋H→TiClx‐（＋HCl（x＝１‐３）のよう
x
１ s）

用されるシャワーヘッドを長寿命化するためのヒントが見に Cl が表面でも脱離することにより Ti 膜が形成されるプ

つかったと同時に，観測や計測が難しい箇所で生じる局所

ロセスを利用していることによる［８］．
従来のプロセスに Ar を加え，更に Ar 流量が H2 流量よ的な現象をシミュレーションによって解明できた意義は非

りも２桁ほど多い Ar ベースのプロセスに関する Denpoh

常に大きい．
ちなみに，CCP のようなタイプの装置であってもガス温らの検討例［８］を以下に紹介する．このプロセスでは，プラ度の上昇を無視できない場合もある［１３］．また，ICP の場

ズマ中の正イオンは大半を Ar＋が占めるようになり，基板合は "!
!比が比較的大きく，一般的にガス温度の上昇を無表面のClの脱離にHよりもAr＋の寄与が大きくなるような視できない［３］．その場合，流体モデルの重粒子のエネル

条件が選ばれたとされる．また，基板に堆積する Ti のス

ギー方程式に，弾性衝突と非弾性衝突によるソース項，おテップカバレッジの改善をめざし，基板に入射する Ar＋のフラックスおよびエネルギーを制御する目的で，シミュレーションにも外部回路を考慮している［８］．以下に Denpoh らの文献［８］で示されている代表的な計算結果を示す．
プラズマ CVD 装置の縦断面を図４に示す．CCP 装置の一つであるが，上部のシャワーヘッドに印加される周波数は 450 kHz と一般的な装置よりも低い．AlN 製のステージ
内には金属のメッシュ電極が埋め込まれており，ステージが下部電極として機能するように設計されている．外部回
路はキャパシタとインダクタを基本とし，下部電極に接続されているメッシュ電極のインピーダンスを制御している．

図３

計算結果の一例を図５に示す．本シミュレーションにも中 心 軸 上 の 電 子 お よ び イ オ ン の 密 度 分 布．（Reprinted
from “Locally enhanced discharges at gas hole outlets
of a showerhead in a plasma etching reactor”, K. Denpoh
2009 J. Phys. D: Appl. Phys 42, 032003 (2009), ! IOP Publishing. Reproduced with permission. All rights reserved）

CFD-ACE＋が用いられており，２４の化学種，８６の気相の反応式，３５の表面反応の式が考慮されている．シミュレー
ションの結果，電極間の主なプリカーサーと正イオンは
TiCl3 と Ar＋となっている．付着係数を比較すると TiCl3
３２９

Journal of Plasma and Fusion Research Vol.93, No.7 July 2017

図４

Ti 製膜を目的としたプラズマ CVD 装置．（Copyright 2016
The Japan Society of Applied Physics）

図６

リアクタンスに対する効果（（a）
上部および下部電極に流れる変位電流 Id と電圧の振幅 Vpp のリアクタンス依存性，（b）
トレンチの底および側壁のステップカバレッジに対するリアクタンス依存性）．（Copyright 2016 The Japan Society of Applied Physics）

が，Ar＋の入射角度に対して強い依存性を持つことに起因している［８］．
図５．
電極間における主な粒子の数密度分布（太い実線は Ar＋）
（Copyright 2016 The Japan Society of Applied Physics）

異なる電力（RF 電源の設定値）において，リアクタンス%に対する Ti 膜のシート抵抗値についても検討されており，投入電力が大きいほどシート抵抗値が下がる傾向が示の方が TiCl2 より小さいと考えられることから，ステップされている．ボーム速度を #% とすると，基板に入射する

カバレッジは大きくなることが期待される［８］．

$Ar＋のフラックスは '!!
(*#% で表される［１６］ことから，電図６に，上部および下部電極で測定された実効的な変位力の増加に伴い Ar＋のフラックスも増加する．また，一般

電流（&'
"
"#
,#
% ），上部電極で測定される電圧の振幅($（$**: Peak-to-Peak Voltage），およびインピーダンスを制に $**も電力に比例するため，電力の増加に伴いイオンエ
ネルギーも増加する．いずれも Ar＋が基板表面に吸着した御する回路のリアクタンス %に対するステップカバレッジ

Cl の脱離を促進していることを示唆している．

を示す．変位電流の式 "!+"$*と $**および "
)$
& !#
(はプラ

以上のように，Ar ベースのプロセスにしてステージにズ マ 密 度 (*と 各 電 極 の シ ー ス 厚 さ )と 関 係 が あ り，

接続されたインピーダンスを制御した結果，Ar＋のフラッ"")となる［１５］．ここで，!+はシース内の電界で*""("(クスが増加して表面の Cl 脱離を促進し，ステップカバレッ

ある．したがって，図６（a）にも示されているように，リア

ジを向上させると同時に，シート抵抗値を低減させて膜質

クタンス %を調整することによりプラズマ密度とシース厚の向上にも成功していることが示された．

さをある程度制御できることがわかる．結果として，図６

FVM を用いた流体モデルでは，化学種や反応式が増え（b）に示されているように，リアクタンス %を調整するこても必要となる物理メモリは一般的なワークステーションとにより，ステップカバレッジを制御できることが示され（workstation ; WS）に搭載できるサイズで十分であり，また．文献［８］では，本シミュレーションとは別のモンテカル

た，計算速度も急激には低下しない．そのため，数多くのロ法を用いたシース内の Ar の軌跡の解析も行い，(*の増化学種と素反応を考慮した計算を安定に解けるのもメリッ

加に伴いシース厚さが減り，イオンと中性粒子との衝突がトの一つである．なお，同じ流体モデルでも有限要素法＋＋

減ることから高エネルギーのAr のフラックスが増加する（Finite Element Method ; FEM）を用いる場合，解が得られことも述べられている（CFD-ACE＋の新しいバージョン

る場合は計算時間を短縮できる傾向がある一方で，問題にでは，壁に入射するイオンのエネルギーおよび角度分布によっては収束しないこともある．また，FEM は一般に膨大

ついて，粒子法によるポスト解析を行う機能が搭載されてな物理メモリを必要とするため，どちらを選択するかは，いるが，文献［８］で述べられているシース内の解析は，これ

解析対象やハードウェアのリソースなども考慮すると良いとは別である）．図６（b）に示されている Bottom Side のスであろう．

テップカバレッジは急峻なピークを持って変化している
３３０Lecture Note

3.2 Processing Plasma Simulation Using Fluid Model３．
２．
４ まとめと今後の課題

K. Ikeda

い．また，表面反応については，ごく限られた条件しかわ

流体モデルを利用したプラズマシミュレーションとし

かっていないのが実情であり，モデルを構築することがでて，Denpoh［４，８］による半導体製造装置のプロセスプラズ

きた場合でも，モデルで必要となる反応速度定数が装置やマシミュレーションを２例紹介した．エッチング装置へのその使用環境に依存するため，汎用的なデータベースはな

応用例では，シャワーヘッドの孔の出口近傍で局所的なプいに等しい．この点は流体モデルに限った話ではないが，ラズマ密度の増加が生じることが解明され，シャワーヘッ

より広範囲で良質なデータベースを準備する必要があり，ドの損傷の原因を明らかにされた．また，プラズマ CVD

今後の大きな課題である．

装置への応用例では，導入ガスに Ar を加え，またその流量

参考文献

を過剰にすることにより，プラズマ中の主たる正イオンが＋Ar となるようにプロセスを変更したことに加え，ステー［１］K. Ono, J. Plasma Fusion Res. 80, 909 (2004).
［２］S.T. Sao et al., Proc. ICCAS 2005, 519 (2005).
［３］K. Ikeda et al., J. Vac. Soc. Jpn. 50, 423 (2007).
［４］K. Denpoh, J. Phys. D: Appl Phys. 42, 032003 (2009).
［５］A. Bhoj et al., ECS Transactions, 25, (8) 719 (2009).
［６］H. Daoxin et al., J. Semicond. 33, 104004 (2012).
［７］M. Miyashita et al., Gasious Electric Conference KW
3.00005 (2015).
［８］K. Denpoh et al., Proc. Syposium on Dry Process, 183 (2016).
［９］S. Hamaguchi, J. Plasma Fusion Res. 85, 177 (2009).
［１０］菅 井 秀 郎：プ ラ ズ マ エ レ ク ト ロ ニ ク ス（オ ー ム 社，
２０００）p.106.
［１
１］菅 井 秀 郎：プ ラ ズ マ エ レ ク ト ロ ニ ク ス（オ ー ム 社，
２０００）p. 48.
［１２］https://en.wikipedia.org/wiki/CFD-ACE+［１３］A. Greig et al., Front. Phys., 2, (80) doi: 10.3389/fphy.2014.
00080.
［１４］O. Takai, J. Plasma Fusion Res. 76, 759 (2000).
［１５］M.A. Lieberman：プ ラ ズ マ／プ ロ セ ス の 原 理 第２版
（丸善，２
０
１
０）p.312.
［１
６］菅 井 秀 郎：プ ラ ズ マ エ レ ク ト ロ ニ ク ス（オ ー ム 社，
２０００），p.57.

ジに接続されたインピーダンスを制御して基板に入射する
Ar＋のフラックスを制御することにより，ステップカバ
レッジと膜質の向上が同時に達成された．どちらもプラズ
マシミュレーションが有効に活用されており，実験のみに頼った従来の方法だけでは成しえなかった例と言える．
FVM を用いた流体モデルによるプロセスシミュレー
ションでは，考慮した化学種や反応式の数が多い場合でも安定して計算を進められる利点がある．また，総セル数が数万程度で収まる二次元の計算モデルであれば，マルチコ
アと 16 GB 程度の物理メモリを有する近年の WS でも十分解析可能である．しかしながら課題もあり，WS とソフトウェアがあれば誰でも簡単にプロセスシミュレーションを行えるかというと，答えは No である．その一番の理由は，シミュレーションで必要とするデータベースの不足が挙げられる．気相の反応，特に複数のガスが混合するプロセス
では，プラズマ中に存在し得る化学種のデータが十分揃っているわけではなく，その素反応についても電子衝突断面積・反応レート等のデータが全て揃っているわけではな!!!!!!!!!!!!!!!!!

だ

けい

池 田

いけ

圭

１
９
８
９年東京理科大学卒業，同年日電アネル
バ㈱（現，キヤノンアネルバ㈱）入社．２
００１
年ウェーブフロント入社．２
０
０
８年１
０月㈱ア
テナシス設立．化学反応を含んだマルチ
フィジックス解析に従事．今回，小田先生より貴重な機会をいただき，また，弊社でサポートしているマルチフィジック
スソフトウェアを応用した代表例として，伝宝様の論文を紹介させていただきました．実用的な応用例として少しでも参考になれば幸いです．

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

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

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

３３１

J. Plasma Fusion Res. Vol.93, No.8 (2017)３
８１‐３８４講座 実践低温プラズマシミュレーション

３．低温プラズマのシミュレーション事例
3. Simulation Examples of Low-Temperature Plasma

３．
３ 流体モデルによる大気圧プラズマシミュレーション
3.3 Atmospheric-Pressure Plasma Simulation Using Fluid Model
小田昭紀
ODA Akinori
千葉工業大学（原稿受付：２
０
１
７年６月１
３日）

先の３．
１節および３．
２節では，低ガス圧力下で生成される低温プラズマの性質を利用した各種材料プロセス技
術のためのプラズマシミュレーション事例に関して紹介された．本３．
３節では，材料プロセスのみならず，医療・バイオ分野への応用展開が精力的にはかられている大気圧下での低温プラズマである大気圧プラズマを対象として，著者の研究グループで開発した低周波駆動されたバリア放電型大気圧プラズマのシミュレーションを例にとり，本プラズマのモデル化，実験との比較を通じた開発したモデルの妥当性検証，本プラズマシミュレーションから得られた成果を順に紹介する．
Keywords:
atmospheric pressure plasma, plasma simulation, fluid model３．
３．
１ はじめに

３．
３．
２ 大気圧グロー放電プラズマとは非平衡大気圧プラズマは，大気圧環境下で熱的に非平衡

大気圧グロー放電プラズマは，上智大学の岡崎らによりな性質（プラズマ中の気体やイオン温度がほぼ室温程度に

１９８０年代後半に初めて見出された日本発のプラズマ源［６］対し，電子温度が数万 K 以上と圧倒的に高い）を有する放である．本放電プラズマの特徴として，（１）減圧環境下で電現象を通じて生成された放電プラズマの総称である．その放電プラズマ生成に必須であった真空装置が必要なくコの際，本プラズマの空間的様相（プラズマ生成が放電空間

ストが低く抑えられる，（２）従来の減圧下で生成された低に対して均一か不均一（局所的）か）や，時間的様相（パ温プラズマと比較して非常に高い化学的反応性を有する，

ルス的か連続的か）などの観点から，大気圧グロー放電，（３）放電電極に対し径方向に均一なプラズマが生成可能なストリーマ放電，誘電体バリア放電，コロナ放電などに大ため大面積に均一な処理が可能である，（４）大気圧環境で別される．これまで，本プラズマは，オゾン生成や産業用

安定なプラズマ生成が可能であるため，反応容器内で放電

紫外光源，材料プロセスへの応用のみならず，現在では農プラズマを生成する従来方式のみならず，大気環境下に放業・水産応用や医療・バイオ分野への応用をめざして研究

出（射出）する方式などプラズマ生成の自由度が高いなど，開発が精力的に展開されてきた［１‐３］．その際，本プラズ

産業応用に優れた特徴を数多く有する．よって，これら特マ技術の発展にプラズマシミュレーションが重要な役割を徴を利用し，材料プロセス（プラスチック材料のぬれ性の果たしてきた．この大気圧プラズマのシミュレーションに向上，液晶基板の洗浄など），医療用器材の低温滅菌，更に関する歴史的経緯や現状等に関しては，当学会誌の小講座

はナノカーボン材料創製など幅広い分野での応用を目指しに詳細が記されているので，そちらを参考にされたい［４］．

研究開発がなされている．本稿では，この大気圧グロー放本稿では，著者の研究グループで開発した低周波駆動さ

電プラズマを低温プラズマシミュレーションの対象としてれたバリア放電型大気圧プラズマのシミュレーション
［５］

扱う．

を例にとり，本プラズマのモデル化，実験との比較を通じ

３．
３．
３ 大気圧プラズマのモデル化た開発したモデルの妥当性検証，その上での本プラズマシ

図１に，著者らのグループで実施された，低周波電圧駆ミュレーションから得られた成果を順に紹介する．

動された誘電体バリア放電型大気圧He/N2 プラズマの計算
モデル図を示す［５］．本図から，両金属電極に誘電体（比誘Chiba Institute of Technology
author’s e-mail: akinori.oda@it-chiba.ac.jp

３８１!2017 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.93, No.8 August 2017

解く必要がある．
!'
"
%#
&
$!$

（３）

!$!%)

（４）

ここで，"，!，)，$はそれぞれ誘電率，電界強度，電位，
空間電荷密度をそれぞれ表す．
本プラズマ中の電子温度は，プラズマ中の原料ガスとの電子衝突の程度を決めるパラメータとして非常に重要である．そこで，この電子温度（もしくは平均電子エネルギー）の時間空間的変化を，以下に示す電子エネルギー保存式から求めている．

図１

低周波駆動・誘電体バリア放電型大気圧 He/N2 プラズマの計算モデル図（文献［５］より引用）．

電率 7.63，二次電子放出係数 0.1）が貼り付けられた誘電体
バリア放電型の電極構成であり，両電極間に形成された

%&
-%"
%'
$%!!%#!!! $,&,!'"
*
$!%#
%
/
,

（５）

#
$%$ ##%!&-%"-*"! "%%&
%'
*'
"
"

（６）

"(
"
%$ ,
%
!

（７）

（放電 プ ラ ズ マ）空 間 に He ガ ス（大 気 圧，300 K，純 度

ここで，"
%，.
%，$,，&,，'"は，平均電子エネルギー，電99.99995％）が封入されているとして，両金属電極間に駆子エネルギーフラックス，各非弾性衝突の閾値，各非弾性

動周波数 15 kHz の低周波高電圧が印加されることで電極衝突の衝突周波数，弾性衝突によるエネルギー損失をそれ間に大気圧プラズマが生成されるとした．その際，不純ガぞれ示す．また，#*，,$，(%は電子の移動度，ボルツマン

スとして N2 ガス（0.5 ppm）が含まれていると仮定した．

定数，電子温度である．

このとき，著者の研究グループで構築した放電プラズマ

大気圧プラズマにおいては，そのガス圧力が非常に高モデルでは，この大気圧プラズマを流体として取り扱い，

く，かつプラズマを構成する電子およびイオン密度が高密

数密度連続の式，電子エネルギー保存式，ポアソン式，熱度になるため，電子およびイオンからの衝突を通じた原料伝導方程式の計４つの支配方程式を空間１次元（電極に対ガスへのエネルギー移行が頻繁に起こる．それにより，プ

して垂直方向）下で数値的に時間空間的に解くことで大気

ラズマ生成用の原料ガスの温度が上昇するため，この影響

圧プラズマの特性を得ている．以下に，これら支配方程式

を適切に考慮する必要がある．そこで，以下の熱伝導方程に関して順に説明する．

式を考慮することにより，大気圧プラズマ中の局所的なガス温度を求めている．

放電プラズマ空間内に存在する電子，正イオン，励起原子・分子の各粒子の時間的振舞いは数密度連続の式から求%($&!' & $%#
&
κ%(&'
!
""#
%
/

まる．数密度連続の式は，放電プラズマ中で考慮した粒子
種 +の密度 -+，粒子フラックスΓ+，粒子が発生する項 #+，消滅する項 %+から構成され，%-+
Γ+"#+!%+
$!%#
%
/

（８）

ここで，$&，!'，(&，κ&，"は，原料ガスの質量密度，定積比熱，ガス温度，原料ガスの熱伝導率，全電流密度をそれぞれ示す．

（１）

数密度連続の式の境界条件に関しては，電子密度およびと表すことができる．ここで，粒子種 +は大気圧プラズマ

各種正イオン密度が両誘電体表面にて表面電荷として蓄積＋中に存在する粒子種（本計算では，e−，He＋，He＋
2 ，N ，するとした．その際，電子に関して，正イオン入射によるN+2 ，He，He＊，He＊
2 ，N，N2 を考慮）を表し，その上でこ放電プラズマ空間への２次電子の放出作用を粒子フラック

れら粒子種が起こす計６１種類の反応過程を考慮した［７‐スの形で考慮した．中性粒子密度に関しては，両誘電体表１０］．このとき，粒子フラックスΓ+は，単位体積・単位時間

面上で全反射する（Neumann 条件）とし，N 原子に関してあたりに通過する粒子数を表し，ドリフト速度 W+およびのみ表面再結合による放電プラズマ空間への N2 の供給作

拡散係数 "+を用いて，用も併せて考慮した．ポアソンの式の境界条件としては，!+$#+-+!"+%-+

接地側金属電極における電位を常時零，駆動側金属電極に（２）おける電位を外部印加電圧とそれぞれ設定した．それに加と与えられる．その際，ドリフト速度が定義されるのは荷え，上述の両誘電体表面への蓄積電荷も両誘電体表面上で電粒子のみであるため，中性粒子に関してはドリフト速度の境界条件として考慮した．
以上の支配方程式を連立して有限体積法を適用して時間が零であることに注意されたい．
電位や電界強度の空間分布を知るためには，空間電荷に空間的に解くことで，荷電粒子や中性粒子密度，電界，電よる電界を考慮する必要した上での以下のポアソンの式を

子温度等の時間空間的変化を求める．その際，粒子フラッ３８２Lecture Note

3.3 Atmospheric-Pressure Plasma Simulation Using Fluid Model

A. Odaクス（（２）式）および電子エネルギーフラックス（（６）式）の数値解法として Scharfetter-Gummel スキーム［１１］を適用した．
上記計算を周期的定常状態に達するまで繰り返すこと
で，大気圧プラズマの時間的空間的な振る舞い，すなわち本プラズマの特性を求めることができる．

３．
３．
４ 大気圧プラズマシミュレーション事例
はじめに，構築した大気圧プラズマの計算モデルが妥当かどうかを検証するために Yuan らの実験（および計算）

図２

低周波駆動・誘電体バリア放電型大気圧 He/N2 プラズマにおける電気的特性の実験および計算との比較；（a）
電流波
形，（b）
電圧−電流特性（文献［５］
および［１
１］
より引用）
．

図３

低周波駆動・誘電体バリア放電型大気圧 He/N2 プラズマ中＋＋の各粒子密度の時間空間分布；（a）
e−，（b）
He2 ，（c）
N2 ，
＋
＊
＊
N，（f）
N2，（g）He ，（h）
He2 （文 献［５］より（d）N4 ，（e）
引用）．

結 果［１２］と 同 条 件（駆 動 周 波 数 15 kHz，印 加 電 圧 振 幅
1.3−2.4 kV）でプラズマシミュレ ー シ ョ ン を 実 施した．
図２に，Yuan らの実験および計算から得られた（a）電流波形および（b）電圧−電流特性を本計算結果と併せてそれぞれ示す．本図（a）から，実験から得られた電流波形に対し，本モデルで得られた波形は，電流ピーク時の値および位相に若干の差異が認められるが，印加電圧の半サイクルに一度の頻度でプラズマの発生および消滅を表す電流パルスの半値幅や，電流パルス消滅後の変化など，本条件下での放電の様子を忠実に再現できている．引き続き，本図（b）において様々な印加電圧（実効値）時における電流ピーク値（実効値）の傾向を調査したところ，本モデルで得られた傾向が Yuan らの行った実験結果と非常によく一致したものとなっていることがわかる．これら結果から，適切なモデル化を行った上でプラズマシミュレーションを実施することで実験結果をよく再現できることが示された．
以上の通りに，構築したモデルの妥当性が検証できた上で，大気圧プラズマの内部構造の詳細を考察していく．
図３に，低周波駆動時の誘電体バリア放電型大気圧 He/N2
プ ラ ズ マ 中 の（a）e−，（b）He＋
（c）N＋
（d）N＋
（e）N，
2 ，
2 ，
4 ，
（f）N2，（g）He＊，（h）He＊
2 の各粒子密度の時間空間分布をそ れ ぞ れ 示 す．こ こ で，本 図 に お け る 0.23 cm お よ び
0.73 cm がそれぞれ駆動側および接地側の誘電体表面の位置を表している．本図（a）の電子（e−）に着目すると，印加電圧の正の立ち上がりにおいて，電極間中央（0.48 cm）と接地側 誘 電 体 表面（0.73 cm）の 間 に お い て 電 子 が最 大2.0×1010 cm−3 もの高密度で急激に生成され，その直後に緩やかに消滅しつつ，印加電圧の変化に応じて駆動側誘電適切に考察することによって，大気圧プラズマ中で起こっ体表面（0.23 cm）から接地側誘電体表面（0.73 cm）に移動ている非常に短時間かつ局所的な現象を詳細に可視化するしている．その後，印加電圧の極性が負で立ち上がった際ことが可能になる．具体的には，両誘電体表面の間につくに，この移動してきた電子を種として，再び電子が急激にられる放電ギャップ間に存在する，大気圧プラズマ中のど高密度で生成され，その後緩やかに消滅している．よって，の粒子種（荷電粒子，ラジカル，励起粒子など）が，どの放電空間内の電子は印加電圧半サイクルに一度毎に極性位置そしてどのタイミング（位相）においてどの程度の密（向き）を変えながらこれら挙動を繰り返し行っているこ

度で生成ならびに消滅を含めた変化をしているかを容易にとがわかる．次に，本図（c）の N＋2 イオン，そして本図（d）

把握することができることになる．

10
−3
の N＋
オー
4 イオンにおいては，放電空間内で常時 10 cm

３．
３．
５ まとめと今後の課題

ダの高密度で存在しており，上述の電子の密度が急激に増加するタイミングに呼応して電極近傍でこれらイオン密度

本節では，大気圧グロー放電プラズマをシミュレーショ

も増加している．その他の図においても，粒子種によってンの対象として，著者の研究グループで実施された低周波

時間的空間的な振る舞い（変化）が大きく異なっているこ

駆動された誘電体バリア放電型大気圧プラズマのシミュとがわかる．

レーションの実践例（モデル化の方法や本プラズマシミュ

これらの結果から，プラズマシミュレーションを実施しレーションから得られる成果など）を紹介した．
３８３

Journal of Plasma and Fusion Research Vol.93, No.8 August 2017

この大気圧プラズマシミュレーションから，適切な数値

間違った計算結果を正しいものと誤解しかねないことに注

モデル化（必要な物理化学現象の考慮），シミュレーショ

意されたい．

ンに使用される妥当な基礎データ（電子断面積，輸送係数，
反応速度定数など），適切な初期条件および境界条件，合参考文献

理的な時間で解を得るための数値解法，これらを考慮に入［１］日本学術振興会プラズマ材料科学第１
５
３委員会：
「大気
圧プラズマ基礎と応用」（オーム社，２００９）．
［２］沖野晃俊（監修）：「大気圧プラズマの技術とプロセス開発」（シーエムシー出版，２０１
１）．
［３］小駒益弘：「大気圧プラズマの生成制御と応用技術」（サイエンス＆テクノロジー，２０
１
２）．
［４］杤久保文嘉 他，プラズマ・核融合学会誌 92, 680 (2012).
［５］小森郷平：千葉工業大学修士学位論文「非平衡大気圧
He プラズマの生成制御に関する研究」
（２
０
１
７）
．
［６］S. Kanazawa et al., J. Phys. D 21, 838 (1988).
［７］X. Yuan and L.L. Raja, IEEE Trans. Plasma Sci. 31, 495(2003).
［８］D. Lee et al., IEEE Trans. Plasma Sci. 33, 949 (2005).
［９］X. Song et al., IEEE Trans. Plasma Sci. 40, 3471 (2012).
［１０］T. Murakami et al., Plasma Sources Sci. Technol. 22,015003 (2013).
［１１］D.L. Scharfetter and H.K. Gummel., IEEE Trans. Electron
Devices 16, 64 (1969).
［１２］X. Yuan et al., Vacuum 80, 1199 (2006).

れてシミュレーションを行うことにより，本プラズマの振る舞い（特性）を定性的のみならず定量的にもおおよそ再現可能なレベルに達していることが示された．
最後に，著者が考える低温プラズマシミュレーションに関する今後の課題を以下に述べる．現状の計算機容量の増
大や CPU の高速化，そして並列計算技術の導入，それらに加え低温プラズマのシミュレータが国内外のソフトウェア
会社から販売されている状況を鑑みるに，プラズマシミュ
レーションを実施するためのユーザー側の障壁は以前よりも低くなっており，またシミュレーションを実施する環境も昔に比べて格段に整っていると思われる．しかしながら，プラズマシミュレーションを自分でプログラムを組ん
で（コーディングして）実施するにせよ，シミュレータを使用して実施するにせよ，低温プラズマに関する知識を持った上でプラズマシミュレーションのモデル化手法やその数値計算技術（数値計算方法，初期条件，境界条件など）などについて，ある程度理解した状況で行われなければ，!!!!!!!!!!!!!!!!!!!!

お

だ

あき

のり

小 田 昭 紀２
００
１年北海道大学大学院工学研究科電子情報工学専攻博士後期課程修了．博士（工学）．現在，千葉工業大学工学部教授．低圧
および大気圧条件下の低温プラズマ（弱電離非平衡プラズマ）の基礎特性解明に関する研究に従事．これまで本学会誌にて数回各種記事を執筆しておりましたが，２
０
１
５年７月から本学会誌の編集委員を拝命しております．当学会誌読者が今回の講座記事を一読いただくことによって，低温プラズマそれ自身，そして低温プラズマシミュ
レーションに関して少しでも興味が湧いてもらえれば非常に嬉しく思います．

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

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

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

３８４

J. Plasma Fusion Res. Vol.93, No.8 (2017)３
８５‐３８６

講座 実践低温プラズマシミュレーション

４．おわりに
4. Conclusion
田中正明
TANAKA Masaaki
ペガサスソフトウェア株式会社（原稿受付：２
０
１
７年６月１
３日）

本講座では，電離度が低く熱的に非平衡である低温プラである，金属表面への窒化やDLC成膜などの表面改質に利ズマのシミュレーションに関し，材料プロセスのための低用される，プラズマイオン注入成膜法のシミュレーション

ガス圧力下の低温プラズマから医療バイオ応用を指向した

では，DC バイポーラパルスを使用した際の生成される各

大気圧下の低温プラズマまで，これらプラズマのシミュ

種炭化水素イオンの入射量などを解析した．その更なる面

レーションの事例を紹介しつつ，この低温プラズマシミュ

内均一性，膜質などの最適化に向け，RF パルスの適用や，レーションを行うことで何がどこまでわかるのかについてパルスの on 時間や duty 比をパラメータとしたシミュレー

説明を行った．以下に，これらの内容をまとめることで，ションが必要であることが示された．３つめの事例である

本講座の締めくくりとさせていただく．

トライボプラズマでは，このプラズマによるトライボケミ

第２章２．
１節および２．
２節では，低温プラズマの定義と熱カル反応のシミュレーションを実施することによって，ト

プラズマおよび高温プラズマとの相違点を述べ，低温プラ

ライボロジーにおける諸問題の解決方法の一助となる成果

ズマが，プラズマ生成部を利用する場合と固体もしくは界が得られた．

面との相互作用を利用する場合とに分け，後者を利用する

本章３．
２節では，流体モデルによるプロセスプラズマシ

各産業分野での例を挙げている．本章２．
３節では，低温プラ

ミュレーションの２つの事例を紹介している．エッチング

ズマシミュレーションの重要性として，特に大気圧プラズ

装置への応用例として，容量性結合型プラズマ装置におけマに関しては，測定の難しさから現象の理解にはシミュるシャワーヘッドの孔の出口近傍のガス圧やプラズマ密度

レーションによる可視化が不可欠であることを述べていをシミュレーションにより求め，計測が困難な箇所での局

る．２．
４節では，粒子モデルによる低温プラズマシミュレー

所的現象が解明され，シャワーヘッドの長寿命化に寄与す

ション手法を述べている．ここでは，ボルツマン方程式のる知見を得ることができた．次に，プラズマ CVD 装置への確率論的解法である，荷電粒子に関する PIC/MC 法，中性

応用例として，TiCl4/Ar/H2 プラズマによる Ti 薄膜生成プ

粒子に関する DSMC 法をそれぞれ解説している．そしてそロセスにおける，インピーダンスおよび電力制御したシ

れぞれの手法に関して長所および短所を述べている．２．
５

ミュレーションを行い，Ar イオンによる基板表面の Cl の節では，流体モデルによる低温プラズマシミュレーション

脱離の促進およびシート抵抗値の低減により，ステップカ

手法として，荷電粒子に関して古くから使用され，実績のバレッジの向上および膜質の向上に寄与する結果が得られあるドリフト‐拡散モデルについて，その特徴を述べていた．
本章３．
３節では，流体モデルによる大気圧下で生成される．また中性粒子に関しては，一流体多成分気体に関するる低温プラズマ（大気圧プラズマ）シミュレーションの事

モデルの特徴を述べている．
第３章３．
１節では，粒子モデルによるプロセスプラズマ例を紹介している．上記プロセスプラズマシミュレーショ

シミュレーションとして以下の３つの事例を紹介していンでの動作圧力と比べ圧倒的に圧力が高いため，それに基る．マグネトロンスパッタ装置のシミュレーション事例でづいたモデル化（高圧力下で顕著に起こる化学反応や，電は，動作ガス圧が 1 Pa 以下と低圧であり，静磁場中での高子と原料ガスとの弾性衝突によるガス温度上昇の考慮な精度な荷電粒子の挙動やターゲットへ入射するイオンの情ど）を適切に行った上でシミュレーションを実施すること報が求められる理由から，粒子モデルにより本装置内におで，定性的のみならず定量的にもよい一致を示す高精度なける荷電粒子および中性粒子の振る舞いの解析が行われシミュレーションが実施できることを示した．

た．その際，スパッタにおいて重要であるターゲット（固最後に，本講座で紹介した低温プラズマシミュレーショ

体）およびイオンとの相互作用の解明が本シミュレーショ

ンを実施するにあたっては，（１）対象となる問題に対してン実施における課題となることが示された．２つめの事例

その物理モデルをよく検討した上で対応する数学モデルを

PEGASUS Software Inc.
author’s e-mail: mtanaka@psinc.co.jp

３８５!2017 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.93, No.8 August 2017

適用すること，（２）その数学モデルを精度良く解くことの座をご一読いただき，その上で必要事項を調べ，低温プラできる数値解法の適用，（３）低温プラズマシミュレーショ

ズマシミュレーションを実施することを通じて，有用な手ンで使用される基礎データ（気相反応の衝突断面積データ

段（ツール）としての低温プラズマシミュレーションの理および速度定数，固体壁での表面反応データなど）を収集

解が更に深まることを願ってやまない．

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

た

なか

まさ

あき

田 中 正 明!!!!!!!!!!!!!!!!

し吟味することが重要である．上記項目を意識しつつ本講!!!!!!!!!!!!!!!!!!!!!!!!!!

１
９
７
４年東京理科大学理学部卒．同年セン
チュリリサーチセンタ㈱（現伊藤忠テクノ
ソリューションズ㈱）に入社，数値解析，
統計解析ソフトウェア開発，核融合炉ダイ

バータプラズマシミュレーションに従事する．２
０
０
２年にペガ
サスソフトウェア㈱を設立し，希薄気体，非平衡プラズマ，
形状シミュレーション等のソフトウェア開発に従事し現在に至る．

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

３８６