Transcript 発表資料 (ppt: 2.9MB)

第19回プラズマエレクトロニクス講習会
「プラズマプロセスの基礎と応用最前線」
‐多様化するプラズマ応用プロセスとその制御‐
プラズマシミュレーション技術
2008年10月30日
ペガサスソフトウェア株式会社
田中正明
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
1
目 次


1. はじめに
2. プラズマシミュレーションに用いられる計算方法











PIC/MCC法
DSMC法
ドリフトー拡散モデル
緩和連続モデル
ハイブリッドモデル
ICP装置における誘導電界
中性粒子シミュレーションとのカップリング
3. 気相シミュレーションの計算例


2.1
2.2
2.3
2.4
2.5
2.6
2.7
3.1 対向ターゲットマグネトロンスパッタ装置
3.2 ドライエッチング装置のシミュレーション
3.3 大気圧沿面放電シミュレーション
4. 形状シミュレーションに用いられる計算方法
5. 形状シミュレーションの計算例


Ar/Cl2プラズマによるPoly-Si エッチングシミュレーション
ボッシュプロセスによるSi エッチングシミュレーション
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
2
1. はじめに
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
3
1.1 プラズマシミュレーションの目的
一般にシミュレーションの目的とは、
●何のために行うのか（Needs）、
●どのような利益があるのか（Benefit)
●シミュレーションはどの程度信頼できるのか（Proof)
プラズシミュレーションの場合、
●プラズマの応用技術として光学的な応用を含むエネルギー分野
●化学的、物理的な応用を含む物質・材料分野
●廃棄物・廃ガス処理などの環境分野
●衛星の推進エンジンへの応用を含む宇宙分野
産業・工業用プラズマはシミュレーション手法（物理モデル）からみたとき、主に
●熱プラズマ
●完全電離プラズマ
●非平衡低温プラズマ
本稿では非平衡低温プラズマを対象としたプラズマシミュレーションに関して述べる。
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
4
1.1.1 何のために行うのか(Needs)
プラズマシミュレーションを行えば何かが判る、という
目的が曖昧な理由からでは、単なる時間の浪費であり、
シミュレーションは意味が無いという結論に陥ることになる。
得られた測定可能な物理量を説明する上で、また
現象のみが分かっているとき、その現象を説明する上で、
何が支配的であるかを知ることは重要である。
従って、シミュレーションにより何が得られれば、その
結果をどのように今後に活かすことができるか、という
問題点と方向性に関するシナリオを持つことにより、
シミュレーションの目的は明確になる
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
5
1.1.2 どのような利益があるのか(Benefits)
● 実験・試作コストの低減を行う
●測定が難しい現象を解明することにより予測の手がかりとなる
●物理量を数値化・可視化することにより現状の把握と仮想実験が可能となる
●運転条件を変化させたときの実験データが得られているとき、
定性的ではあっても、目的とするものの最適化が可能となる。
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
6
1.1.3 どの程度信頼できるのか(Proof)
●プラズマシミュレーションには気相中、構造物表面での物理的、化学的
そして光学的な実現象を正確に取り入れるモデルは確かに不可能に近い。
なぜなら、気相中での複雑な反応に伴う基礎データの欠如、もしくは信頼性
そして取り扱う境界条件が時々刻々変化している場合が多いからである。
しかし、
●運転条件、装置条件に適合した、理論的な背景に基づくモデル
によるシミュレーションは、実現象を模擬したものであり、
得られる数値も定性的には信頼に値するものである。
●物理現象が比較的明確な対象によっては定量的に信頼できる
結果が得られる。
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
7
2. プラズマシミュレーションに用いられる計算手法
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
8
大別して
●気相中でのプラズマ（荷電粒子）とラジカル・励起種・供給ガス（中性粒子）の
挙動をシミュレーションする計算手法
●基板における微細形状の時間変化をシミュレーションする計算手法
プラズマ装置内の気相シミュレーションにおける解析領域と
微細形状のシミュレーションにおける解析体領域の比は10５以上になる。
近年常圧下での数mm以下の気相シミュレーションも頻繁に行われる
ようになってきた。
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
9
気相シミュレーションにおける計算手法の種類
電磁場の基礎方程式：マクスウェル方程式 と
気体運動論の基礎式：ボルツマン方程式 を
どう解くか？
●粒子モデル
●流体モデル
●ハイブリッドモデル（粒子モデルと流体モデルの混合型）
●グローバルモデル
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
10
粒 子 モ デ ル
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
11
2.1 PIC/MCC法(荷電粒子のための粒子モデル）
PIC/MCC法は、Particle-in-Cell(PIC)法と Monte Carlo Collision(MCC)法を
組み合わせてプラズマ（荷電粒子）の挙動をシミュレーションする手法である。
PIC（Particle-In-Cell）法はおよそ40年以上前に流体方程式の解法として開発。
その後核融合プラズマにおける電磁流体方程式解法にも使用され、
およそ20年前から弱電離プラズマ（非平衡低温プラズマ）シミュレーション
に使用され、今日に至っている。
計算手法の詳細は例えば Birdsall の教科書[1]、南部[2]などを参照して
頂きたい。
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
12
計算手順の概要は次のようなものである。
1. 静電位に関するポアソン方程式を解く
ポテンシャル計算のために各セルの電荷密度を求め、電位の境界条件を
考慮した上で各セルのポテンシャルと電界を計算する。
2. 荷電粒子の運動方程式を時間積分する
セル境界上の電界、 磁界を粒子位置に内挿して Δt 後の粒子の速度と
位置を計算する。
3. 壁面での境界条件
壁面での荷電粒子の消滅や散乱、電荷の蓄積、2次電子放出など
境界上での粒子の計算を行う。
4. 中性粒子との衝突
荷電粒子と中性原子・分子間の衝突の種類 (弾性散乱、電離、励起、解離)
を衝突断面積に従った確率で決定し、衝突後の速度を計算し、電離による
イオンと電子の生成も考慮する。
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
13
この手順を、プラズマの密度分布や電位分布が(周期)定常状態になるまで
短い時間ステップ Δt で繰り返していく。
前記の 1～3 を繰り返して計算するのが PIC法であり、4 が MCC 法である。
PIC/MCC法では、例えば 108 ～ 1010 個の荷電粒子をまとめて、1個の超粒子
として扱い、この超粒子の運動を追跡する。
メッシュ分割した解析空間内に電子とイオンの超粒子を配置し、その超粒子の
分布から決まる空間電荷をグリッド上に振り分ける。そして与えられた境界条
件(電極の電圧やアースの配置など)の下でポアソン方程式
を解いて電位Φを計算し、その勾配から電界を求める。
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
14
荷電粒子の運動方程式は、
ここで、 m: 荷電粒子の質量; E: 電界;B: 磁界; q: 荷電粒子の電荷;
r: 荷電粒子の位置座標; v: 荷電粒子の速度
この運動方程式を、良く知られている蛙飛び法、ベルレ法などで定式化し
Δt 後の位置と速度を計算する。
PIC/MCC法における留意点
1）空間電荷の平滑化
2）ポアソン方程式解の数値誤差
3）メッシュ幅
4）超粒子の重み
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
15
2.2 DSMC法(中性粒子のための粒子モデル）
プラズマ装置内の原料ガスの流れ場などのように圧力が低い流れ場について
解析する場合、もしくは常圧下であっても解析体系がμmとなるとき、気体を
連続体であるとみなして取り扱うナビエ-ストークス方程式は適用できず、
Boltzmann方程式にまでさかのぼって考えなくてはならない。
ところが、この Bolzmann方程式は複雑な非線形積分微分方程式であり、
その取扱いは非常に難しい。
DSMC法 (Direct SimulationMonte Carlo method、 直接シミュレーションモンテ
カルロ法)はBoltzmann方程式を直接解くのではなく、Boltzmann方程式の
基礎になっているそれぞれの粒子の衝突過程を確率的に取り扱うことによって
流れ場を解析するといった方法である。
この手法の詳細は、南部[2]、Bird [3] そして南部[4]の教科書を参照して頂きたい。
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
16
DSMC法では、気体分子をPIC/MCC法と同様、超粒子として取り扱う。
真空装置内といえども装置内の全ての気体分子を計算機で追跡すること
は不可能で、例えば 10１５ ～ 1020 個の気体分子をまとめて1個の超粒子
として扱う。
数万～数100万個程度の超粒子を取り扱い、各粒子の位置、速度、内部状態等
をメモリー内に記録し、衝突や境界の影響によってそれらの値を更新していく
という手続きを繰り返す。
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
17
分離の原理
●メッシュ幅が平均自由行程程度
●クーラン条件
が満たされているとき、粒子間衝突と粒子の移動が分離できる。
衝突過程については、与えられた衝突断面積に従って確率的に引き起こし、
衝突後の2粒子の速度もその衝突物理モデルに従い確率的に決定する。
具体的には、
流れ場を適当な大きさのセルに分割し、 同一セル内の 2 個の粒子を物理モデ
ルに基づいた確率則に従って衝突対として選択するというものである。
（現在では、最大衝突数法が最良）
衝突後の粒子対の速度は、エネルギー、運動量等が保存するように決める。
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
18
DSMC法により直接得られるものは、多数の粒子の位置、速度等のサンプル
データである。
このデータを加工→粒子の速度分布関数、密度分布、流速分布、温度分布
統計誤差→サンプル数の1/2乗に反比例
時間平均のゆらぎが小→十分な精度を持つまでサンプリング
プラズマ装置内もしくは真空装置内では、密度差が大きく異なる
原子・分子が存在し、小さな密度を持つ原子・分子の挙動を正確にシミュレー
ションするには、従来のDSMC法ではすべての粒子種の重みは同一であるた
め、膨大なサンプル数が必要であるが、Weight algorithm(南部[2])を採用
することにより、ほぼ同数のサンプル数で行える。
DSMC法の留意点
1）メッシュ幅
2）サンプリング
3）分子模型（剛体球、VHS、VSS）
DSMC法のトピックス
・圧縮性流体解析のための新手法
（メッシュ幅>>平均自由行程）
粘性流～自由分子流
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
19
流 体 モ デ ル
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
20
2.3 ドリフト－拡散モデル（荷電粒子のための流体モデル)
ガス圧が高くなると荷電粒子が分子と衝突する際の平均自由行程が
小さくなる→PIC/MCC 法（粒子モデル）を用いようとすれば→
空間のメッシュ分割幅と時間進展の Δt をともに小さくする
必要があり  計算時間が実用的な範囲を超えてしまう。
そのような場合には、電子やイオンの集団をそれぞれ連続体(流体）として
取り扱うモデルを用いる。
流体モデルでは、電子とイオンの各々についてマクロな量としての
数密度
、流速ベクトル
(r は位置座標、t は時刻）
を未知変数として記述される
連続の式
運動方程式
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
21
を電位のポアソン方程式とカップルして解くことにより、プラズマ密度や
ポテンシャル、および電子・イオンのフラックスを求める。ここで
添え字の α は荷電粒子の種類を識別する指標で、 α=e は電子を、
α=ion はイオンを意味する。
また必要な場合には電子の平均エネルギー についての
バランス方程式
を解く。
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
22
基礎方程式の導出
基礎方程式（粒子種を表す添字は省略）
∫
1
c
eZ
∂f
+ c・▽rf +
E・ ▽cf
m
∂t
c2/2
=
1
c
∫
c2/2
∂f
∂t
( )
col
d 3c
d3c
E ：電界
∂f
∂t
( )
col
：衝突項
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
23
基礎方程式
連続の式
∂n
+ ▽・nv = R
∂t
運動方程式
∂mnv
∂t
+ ▽・mnvv + ▽・P - neZE= R
mom
エネルギーバランスの式
∂n<ε>
∂t
+ ▽・nv<ε> + ▽・P ・v + ▽・Q + env・E= R
eng
R:電離等の反応による粒子の生成レート
Rmom:衝突による運動量の交換レート
Reng:衝突によるエネルギーの交換レート
P:応力テンソル=mn<(c-v)(c-v)>=mn<c’c’>
Q:熱流ベクトル=mn<c’2c’>/2
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
24
生成レート（電子の場合）
R= (νi – νa) n
Rmom= - mnνm v
Reng= - nνm(2m/M) <ε>
- nνiEi –nνexEex - ……
νk/N = ∫σk(c)f(c)c d3c
νi
νa
νm
νex
Ei
Eex
N
M
σk
：電離の衝突周波数
：電子付着の衝突周波数
：運動量移行衝突周波数
：励起の衝突周波数
：電離エネルギー
：励起エネルギー
：衝突相手分子の数密度
：衝突相手分子の質量
：k衝突の衝突断面積
（分子は電子に対して静止しているという近似）
・各種衝突について、衝突断面積のデータが必要
・電子の速度分布関数を知る必要
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
25
運動方程式の簡略化（ドリフトー拡散モデル）
・応力テンソルPは等方的  スカラー圧力 p
p=mn<c’2>/3 = 2nεT/3 = nkBT
1
mnνm
T:電子温度
∂mnv
+ ▽・mnvv + ▽nkBT + neE = - v
∂t
左辺第１項のオーダー : v/νmτv=(τm/τv)v
左辺第２項のオーダー: v2/νmL = (τm/L)v2 = (λm/L)(v/vT) v
τm :運動量緩和に関する平均衝突時間（1/νm）
τv : マクロな速度が変化する時間尺度
λm : 運動量緩和に関する平均自由行程
L :マクロな速度が変化する空間尺度（例えば電極間距離）
vT: 電子の熱（ランダム）運動の速度スケール √<c’2>
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
26
運動方程式の簡略化（ドリフトー拡散モデル）
・τm/τv ≪ 1 の場合は第１項（時間微分）は無視できる。
アルゴン中の電子の場合 νm ～ 1010ptorr [/s]
τv=107 [s]（10MHzに相当）とすると
τm/τv～10-3/ptorr
p > 10 [mTorr] なら条件が成り立つ。
・v/vT ≪ 1 の場合は第２項（対流項）は無視できる。
プロセスプラズマではだいたいこの条件は
満たされていると考えてよい。
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
27
運動方程式の簡略化（ドリフトー拡散モデル）
圧力勾配項
1
nevT2
1
▽nekBTe
meneνm = ne ▽( 3νm ) = ne ▽(
1
De
= ne ▽(neDe)～
ne ▽ne
eE
meνm
外力(電界）項
neτmvT2
)
3
= μeE
μe =
e
meνm
拡散係数
De =
τmvT2
3
フラックス
Γe = - neμe E - De▽ne
移動度
=
kBTe
meνm
=
μekBTe
e
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
28
エネルギーバランス（電子）
∂n<ε>
∂t
+ ▽・nv<ε> + ▽・P ・v + ▽・Q + env・E= R
eng
・熱流ベクトルは温度勾配に比例
Q= - K▽Te = -
2K
▽εT = 3kB
5
neDe▽εT
3
(K=5De/2)
・圧力はスカラー
▽・P・ve = ▽・(pve) = 2 ▽・(neεTve) =
3
2
▽・(ΓeεT)
3
・平均速度の運動エネルギーは、熱運動のエネルギーに対して無視
<ε> ≒ εT ≫ meve2/2
・エネルギー損失
Reng = -ne(νm
2me ν
+ i
M
Ei
Eex
)<ε> = -neνε<ε>
+νex
<ε>
<ε>
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
29
電子に関する基礎方程式のセット
連続の式
∂ne
∂t
運動方程式
Γe = - neμe E - De▽ne
+ ▽・Γe = (νi -νa) ne
エネルギーバランス式
∂
∂t
5
(ne<ε>) +▽・ 3 Γe<ε>
-
5
neDe▽<ε>
3
+eΓe・E
= -neνε<ε>
・νε≪ νm であるのでエネルギーバランス式においては
一般に時間微分項を無視することはできない。
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
30
イオンに関する基礎方程式のセット
・運動方程式に関しては、電子と同様ドリフトー拡散モデル
・イオンのエネルギーは電離や電子付着にほとんど影響しない
ので、エネルギーバランス式は解かない。イオン温度は一定
と仮定。
連続の式
∂nion
+ ▽・Γion =
∂t
νine ( 正イオン)
νane ( 負イオン)
運動方程式
Γion = ± nionμion E - Dion▽nion
μion =
(mion+M)e|Z|
mionMνm,ion
Dion =
μionkBTion
e
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
31
計算手順
これらの基礎方程式を離散化して解く
E
電位に関するポアソン方程式
Γe、Γion を計算
連続の式
ne, nion
電子エネルギーバランス式
t=t+Δt
<εe>
この間、中性粒子種の
密度は変化しないと考
える
ICPのとき、
コイル電流から
生ずる誘導電界
・パワー
モンテカルロ法
による電子エネ
ルギー分布関数
の計算
ハイブリッドモ
デルの場合
中性粒子種に
関する流れの
計算
数周期毎に１回やり取りする
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
32
電子の速度分布関数に対する仮定
衝突周波数： νk/N = ∫σk(c)f(c)c d3c
（１） 局所マクスウェル分布を仮定
f(c)=fLM(c;Te)= 4π
me
2πkBTe
3
2
2
m
c
e
c2 exp 2kBTe
νk/N はTeの関数
色々なTeの値に対してνk /N を計算し、Te～ νk /Nの関係を
テーブルとして用意しておく。
連続の式
運動方程式（粒子フラックス）
電子エネルギーバランス式
電位のポアソン方程式
Te（セル毎）
νk
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
μe, De
R, Rmom,Reng
33
電子の速度分布関数に対する仮定
（２） 一様電界中の電子スウォームの分布関数を用いる
一様な電界Eの中を移動する電子群に対する２項近似の
ボルツマン方程式を解き分布関数を求める。
Nμe、NDe、νk/N および<ε>はそれぞれ換算電界E/Nの関数
色々なE/Nの値に対して計算を行い、 Nμe、NDe、νk/N ～E/N
またはNμe、NDe、νk/N ～ <ε>関係のテーブルを作成しておく。
連続の式
運動方程式（粒子フラックス）
電位のポアソン方程式
E/N（セル毎）
μe, De
E/Nの関数
R, Rmom,Reng (局所電界近似)
または
電子エネルギーバランス式
<ε>（セル毎）
μe, De
<ε>の関数
R, Rmom,Reng
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
34
イオンの輸送パラメータ
電子と異なる点
・質量 mion ～ M
・分子速度が無視できない
νm/N = ∫∫σm(cr)f(cr)cr fM(cn)d3c d3cn
cr=c-cn , cnは分子の速度
・分子は多様なので、衝突断面積のデータがみつから
ない場合も多い
ドリフト速度 vd は換算電界E/N に依存する。
vd ∝ E/N (E/Nの小さい領域 ≪ 5.6 Td )
(1Td=10-17 volt・cm2 )
∝√(E/Nの大きい領域) (E/N  ∞)
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
理論が難しい
35
イオンの輸送パラメータ
（１） 実測データを使用する
・vd, Nμ,ND ～E/N の関係をテーブル表にしておく。
・文献等ではμそのものではなく、μ0 (reduced
mobility )の値をE/N またはE/p0に対してプロット
していることが多い。
μ0=(p0/760)μ, E/N[Td] =2.828×E/p0[volt/cm/Torr]
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
36
イオンの輸送パラメータ
（２） 次の近似式からμを求める。
3e|Z| π(mion+M)
Nμ=
8
2mionMkBTeff
1/2
1
Ω(1,1)(Teff)
kBTeff=kBTn +(1/3)Mvd2
Ω(1,1)(Teff)=
1 (k T )-3∫ε∞2σ (ε)exp(-ε/k T )dε
B eff
m
B eff
0
2
D=
μkBTeff
e
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
37
境界条件
・境界面における粒子フラックスを規定する。
Γdrift =
(ξnμE・n)^n ^（ξE・n >^0 のとき）
0
（ξE・n^> 0 のとき）
Γdiff = (-D ▽n・n^) n^
（n境界=0）
Γout = aΓdrift + bΓdiff
(a,b=0 または1)
ξ=±1
あるいは
Γout = (1/4)nvT n^, vT=√(8kBT/πm)
n: 境界面における外向き単位法線ベクトル
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
38
境界条件
通常、電子に関しては
Γe,out = (1/4)neveT n^
Γe,in = - Σαγse,αΓα,out
γseは正イオンの衝突による２次電子放出係数。
総和は正イオン種αに関して取る。
Γe,境界= Γe,out + Γe,in
イオンに関しては
Γion,境界= Γion,out = Γion,drift + Γion,diff
電子エネルギーフラックスに関しては
qe,境界=2kBTeΓe,out + 2kBTe,in Γe,in （Te,inは仮定）
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
39
時間刻みΔt
Δt < min Δr/(μeEr), Δz/(μeEz)
・電界ドリフトによる移動が１ステップで１メッシュを
こえない。
・電子のドリフト速度で評価
・生成レートが大きい場合には、この評価よりもかなり
小さくしなければならない場合もある。
・およそΔt =10-10～10-12 程度
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
40
2.4 緩和連続（RCT)モデル（荷電粒子のための流体モデル)
・ 連続の式の粒子フラックスにドリフト速度を用いる。
・ 運動量保存則と電子のエネルギー保存則に
それぞれの緩和時間を取り入れた真壁等によるモデルがある。
(文献[6],[7],[9])
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
41
2.5 ハイブリッドモデル（流体モデルと粒子モデルの混合型)
1） 電子のモンテカルロシミュレーションを併用する
電子の速度分布関数（エネルギー分布関数）についてのみ PIC/MCC 同様の
手法を用いて計算セルごとに値を求める。
ただし、この場合にはプラズマ密度の空間分布によってプラズマ中の電界が
変化するので、流体モデルによるプラズマ密度および電界の計算と
モンテカルロ計算を数周期 ～ 数十周期ぶんずつ交互に行いながら計算。
計算時間が多少増えるが、分布関数を仮定する必要が無いのが利点である。
電子エネルギー分布関数から各セルにおける電子温度、電子衝突による
反応レートが求まる。
（文献[8]）
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
42
プラズマの流体モデルの計算で得られた電子密度、分子密度、
および（一般には時間依存の）電界E(r,t)をもちいて、
同じ計算体系（形状）で電子のモンテカルロシミュレーション
をおこなう。
電子の速度分布関数f(r,c)をえる。（セル毎）
（ 実際にはエネルギー分布関数f(r,ε) ）
νk および Te がセル毎に求まる。
連続の式
運動方程式（粒子フラックス）
μe, De
R, Rmom,Reng
ne(r)（周期平均）
電界E(r,t)
電位のポアソン方程式
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
43
2） 電子をPIC/MCCで計算しイオンは流体モデルで計算する
マグネトロンスパッタ装置におけるマグネトロンプラズマシミュレーションは
PIC/MCC法のみで計算する方法と電子のみPIC/MCC法を用い、
イオンの生成率を基にしたイオンに関する流体モデルによりイオン密度を
求める方法がある。マグネトロンプラズマに限らず、電子の非線形性が強く、
流体モデルで近似できない場合に採用される。
（文献[9]-2)
[9]-2 志堂寺 栄治 、真壁 利明
マグネトロンプラズマのシミュレーション
真空, Vol.43, No.8, pp. 767-772, 2000.
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
44
3） 電子を流体モデルで近似しイオンの挙動はモンテカルロ法を用いる
ECRプラズマシミュレーションにおいて磁力線に沿って電子の流体モデルの
定常解を求め、イオンの密度はモンテカルロシミュレーションにより求める方法。
（文献[10]）
4）シース領域のみ電子のモンテカルロシミュレーションを併用する
2次電子放出によりプラズマが生成・維持されるとき、シース領域における
電子は強い非平衡性と高エネルギーを持つため、モンテカルロ法による
シース端での生成率を、バルクプラズマにおける流体モデルへ引き渡す。
（文献[9]の他、1990代初頭から多くの論文）
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
45
2.6 ICP装置における誘導電界
ICP（軸対称モデル）の場合
▽2Eθ – iωμ0σEθ= 0 を解く
Eθ：誘導電界
ω：コイル電流の角振動数
μ0:真空の透磁率
σ：プラズマの導電率
連続の式
運動方程式
=
e2ne
me(νm+iω)
ne(r)（周期平均）
νm
電位のポアソン方程式
Eθ
σ|Eθ|2(Power)
電子モンテカルロ計算
（ハイブリッドの場合）
電子エネルギーバランス式
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
46
ファラデーの法則
アンペールの法則
ガウスの法則
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
47
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
48
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
49
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
50
● イオンモンテカルロシミュレーション
今まで述べた流体モデルではイオンエネルギーを一定として解いている。
しかし、一般にシース領域ではイオンは加速され高エネルギーとなる。
従ってシースを通過して基板、ターゲットへ入射するイオンの速度分布
関数を求めるために、流体モデルもしくはハイブリッドモデルから得られた
プラズマパラメータおよび時空間電界を用いてシース内のイオンの挙動を
モンテカルロ法により計算し、入射角度・エネルギー分布関数を求める。
●グローバルモデル
一般にゼロ次元のレート方程式と電子エネルギーバランス式を用いて
時間依存の初期値問題としてスティフな系での常微分方程式を解く。
長所は、イオン種、ラジカル種が多く存在存在するとき、短時間でおよそ
の絶対値が把握出来る。
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
51
2.7 中性ガスシミュレーションとのカップリング
衝突による励起、解離
ガスの流動
原子、励起種生成
原料ガスの消費
中性粒子種の密度変化
電離、電子付着などに影響
・流体モデル
各種成分の質量保存式（対流－拡散形） nα （α:中性粒子種）
混合ガスの運動方程式
vgas
混合ガスのエネルギー保存式 Tgas
・希薄気体（DSMC法）
・電子の時間スケールとガス流れの時間スケールは
O(104)以上異なるので分離し、交互に計算を行う。
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
52
3. 気相シミュレーションの計算例
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
53
3.1 対向ターゲットマグネトロンスパッタ装置
粒子モデル
PIC/MCC
DSMC
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
54
解析モデル
残留磁束密度：1 T
N
Ar,O2:
0.25 Pa
200
S
-400 V
S
N
80
100
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
55
磁力線と磁束密度
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
56
荷電粒子密度
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
57
O(-)密度、電子、電子温度
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
58
Ar(+),O2(+)入射フラックス、Siスパッタ粒子フラックス
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
59
Si 密度、温度
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
60
Si 速度、堆積フラックス
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
61
3.2 ドライエッチング装置のシミュレーション
参考文献：M. Shiozawa and K. Nanbu,
"Coupling of plasma and flow in materials processing",
Thin Solid Films, Vol.457,pp.48-54,2004
ハイブリッドモデル
＋
ＤＳＭＣ
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
62
解析モデル
121
100
RF
Ｃｌ２
ICPパワー：500[W]
13.56[MHz]
Ｃｌ２流量 ：50[sccm]
385
z
r
250
15
78
174
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
63
反応式
＜電子衝突断面積＞
e-+Cl : 弾性、非弾性（励起、電離）等３種
e-+Cl2 : 弾性、非弾性（励起、電離、解離、解離性付着）等６種
＜反応レート定数＞
Cl- + Cl2+ → ３Cl
Cl- + Cl+ → Cl２
Cl+ + Cl2 → Cl2+ + Cl
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
64
コイルパワー吸収分布、誘導電界強度分布、
電子エネルギー分布関数
EEDF
（電子エネルギー
分布関数）
×
×
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
65
電子に関する物理量、電位分布
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
66
Cl-、Cl+密度分布、 Cl-、Cl+生成率分布
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
67
Cl2+密度、生成率分布、 Cl2、 Cl分布
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
68
3.3 大気圧沿面放電シミュレーション
流体モデル
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
69
解析対象
図のような2次元体系で、電極間距離 d (0.5mm ～5mm)、
誘電体厚み h (0.3mm,2mm)と放電開始電圧V の
関係を調べる。
仮想境界
仮想境界
窒素 1 [atm], 300 [K]
陰極
h
1
d
1
0.3
plasma
仮想境界
0.3
陽極
誘電体 比誘電率εr=3.0
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
印加電圧
V[V]
70
解析モデル形状
N2 1[atm], 300[K]
1
d
1
1
陽極：電圧V
陰極（接地）
誘電体 εr=3
0.3
h
単位：[mm]
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
71
ガス条件
・ガス条件 ： 窒素
圧力 1 [atm]、 温度 300 [K]
・粒子種 ： N2、N、N2*（励起種）、e-、N2+
・考慮した衝突反応 ：
e- + N2  e- + N2 （弾性衝突）
e- + N2  2e- + N2+ （電離）
e- + N2  e- + 2N （解離）
e- + N2  e- + N2* （振動励起）
e- + N2  e- + N2* （回転励起）
ガス圧および温度は一様かつ一定と仮定する。
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
72
境界条件
・静電ポテンシャルに関する境界条件
[陰極] φ= 0 [V]
[陽極] φ= V [V] （V の値を変化）
[誘電体の下面] φ= 0 [V] （一部遠方接地）
[仮想境界] 当該境界面に垂直な方向に 100 [mm]
離れた場所でφ= 0 [V] に規定（遠方接地）
・電子フラックスに関する境界条件
Γe ・ n= (1/4) vT ne - γ Γi・n
n は外向き単位法線ベクトル、vT = √(8eTe/πme)、
γはイオン衝撃による２次電子放出係数。
誘電体および電極表面ではγ=0.05、その他の境界ではγ=0
・イオンフラックスに関する境界条件
Γi ・n= (ni μi E – Di ▽ni )・ n ただし境界において ni = 0
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
73
初期密度と輸送パラメータ
・電子/イオンの初期密度: ne = ni = 108 [/m3]
・輸送パラメータ （移動度、 拡散係数、レート係数）
[電子] : あらかじめ定常タウンゼント放電に対する２項近似の
ボルツマン方程式を解いて、移動度 を換算電界 E/N
の関数として、また拡散係数および各種電子衝突の
レート係数は電子温度（平均エネルギー）の関数として
テーブルを作成し、それをもちいて補間する。
[N2+] : 文献[5] の実測データから移動度 および拡散係数を
E/N の関数としてテーブルを作成し、それをもちいて
補間する。
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
74
シミュレーション結果 (h=2mm)
放電開始電圧－電極間距離 (h=2mm)
4.0
3.5
電圧 [kV]
3.0
2.5
2.0
1.5
2mm
1.0
0.5
0.0
0.0
1.0
2.0
3.0
4.0
5.0
6.0
電極間距離 [mm]
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
75
シミュレーション結果 (h=0.3mm)
放電開始電圧－電極間距離 (h=0.3mm)
1.6
1.4
電圧 [kV]
1.2
1.0
0.8
0.6
0.3mm
0.4
0.2
0.0
0.0
1.0
2.0
3.0
4.0
5.0
6.0
電極間距離 [mm]
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
76
シミュレーション結果 (遠方接地)
放電開始電圧－電極間距離（遠方接地）
7.0
6.0
誘電体表面において
表面電荷の蓄積を考
慮せず
電圧 [kV]
5.0
4.0
3.0
2.0
0.3mm
1.0
100mm
0.0
0.0
1.0
2.0
3.0
4.0
5.0
6.0
電極間距離 [mm]
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
77
シミュレーション結果 （ｈ依存の比較）
放電開始電圧－電極間距離
14.0
12.0
■ パッシェン則（空気）
電圧 [kV]
10.0
■ 遠方接地
h2
◆mis
h=2
h0.3
▲Pas
h = 0.3
8.0
6.0
4.0
2.0
0.0
0.0
1.0
2.0
3.0
4.0
5.0
6.0
電極間距離 [mm]
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
78
プラズマ成長／消滅の判断（h=2mm, d=4mmの例)
密度ピーク値の
時間変化
放電
消滅
電子密度分布
V=3200[V]
V=3300[V]
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
79
密度ピーク値の経時変化（h=2mm）
d=0.5mm
V=1200V
d=3.0mm
V=3200V
d=1.0mm
V=2100V
d=4.0mm
V=3300V
d=1.5mm
V=2700V
d=2.0mm
V=3100V
d=5.0mm
V=3500V
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
80
電子密度の空間分布（h=0.3mm）
d=1.0mm
V=1500V
d=4.0mm
V=1400V
d=2.0mm
V=1500V
d=3.0mm
V=1400V
d=5.0mm
V=1500V
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
81
電子密度の空間分布（h=2mm）
d=1.0mm
V=2100V
d=4.0mm
V=3300V
d=2.0mm
V=3100V
d=3.0mm
V=3200V
d=5.0mm
V=3500V
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
82
4. 形状ミュレーションに用いられる計算手法
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
83
表面反応モデルと形状進展モデル
表面反応モデル：
粒子の輸送、堆積、エッチング、スパッタリング、
そして反射に関するモデル
・エッチングレート、デポジションレートを求める
・素反応をモンテカルロ法などで計算
形状進展モデル：
・ストリングモデル（格子点の移動）
・レベルセット法（表面曲線の時間発展方程式を解く）
（慶大：真壁研究室）
・セル法（セルの生成、消滅）
（京大：斧研究室、ミシガン大：Kushner等のグループ）
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
84
モンテカルロ法における反射確率と反応確率
鏡面反射確率（f） = (E/ Eth )(θ- θcutoff)/(90- θcutoff)
ここで、 E、 θは入射エネルギー[eV]、法線となす入射角度[度]、
そして、鏡面反射エネルギー＝E*f として計算。
fd を拡散反射確率として、U ： f / (f + fd ) で決定
また、表面反応が生じないとき、反射モデルとして使用。
反応確率の種類は、
1:
2:
3:
4:
5:
P0
P0 (E/E0 –1) 1/2
P0 max(0, 1 - E / Ec )
P0 (En - Ethn ) / (Ern - Ethn ) f(θ)
p0 (E0.5 – Eth0.5 ) f(θ)
（パラメータ： P0 ）
（パラメータ： P0,, E0 ）
（パラメータ： P0 , Ec ）
（パラメータ： P0 , Eth, Er,, n, f(θ) ）
（パラメータ： p0 , Eth , f(θ) ）
E, Eth, E0, Ec, そしてErの単位は[eV]、 f(θ) は入射角度依存の相対確率
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
85
5. 形状シミュレーションの計算例
表面反応モデル：素反応をモンテカルロ法
形状進展モデル：セル法(固体層占有率、表面被覆率)
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
86
Ar/Cl2プラズマによる
Poly-Si エッチング
参考文献：
R. J. Hoekstra, M. J. Grapperhaus and M. J. Kushner,
"An Integrated Plasma Equipment Model for Polysilicon Etch Profiles in
an Inductively Coupled Plasma Reactor with Subwafer and Super wafer Topography“
J. Vac. Sci. Technol. A. 15, 1913 (1997).
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
87
Poly-Si：Ar/Cl2プラズマ表面反応機構(１)
Cl2+
Cl+
SiCl2
Cl2+
Si
Cl
SiCl2
SiCl2
SiCl
Cl+
Cl2+
バルク表面種
SiCl2
錯体
Cl SiCl2
Cl2+
SiCl2
SiCl2
SiCl2
SiCl2
SiCl2
Ar+
SiCl3 ,Ar
Ar+
SiCl4
デポ反応
イオンアシスト
SiCl2 ,Ar
エッチング反応
Cl
I*
SiCl3
Cl,Cl+
高速中性粒子
SiCl4
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
88
Poly-Si：Ar/Cl2プラズマ表面反応機構(2)
p0
Si(s) +Cｌ(g)  SiCｌ (s)
SiCｌ(s) +Cｌ(g)  SiCｌ2 (s)
SiCｌ2(s)+Cｌ(g)  SiCｌ3 (s)
SiCｌ3(s)+Cｌ(g)  SiCｌ4 (g)
Si(s) +SiCｌ2(g)  Si(s)+ SiCｌ2 (s)
SiCｌ(s) +SiCｌ2(g)  SiCｌ(s)+ SiCｌ2 (s)
SiCｌ2(s)+SiCｌ2(g)  SiCｌ2(s)+ SiCｌ2 (s)
SiCｌ3(s)+SiCｌ2(g)  SiCｌ3(s)+ SiCｌ2 (s)
0.99
0.20
0.15
0.0001
0.8
0.5
0.3
0.1
1/2
SiCｌ2(s)+Ar+(g)  SiCｌ2(g)+ Ar (g)
0.16 ( E/E0 -1 )
SiCｌ3(s)+Ar+(g)  SiCｌ3(g)+ Ar (g)
1/2
E/E
0.16
0( -1
)
1/2
0.13 ( E/E0 -1 )
SiCｌ(s) +Cl+(g)  SiCｌ2(g)
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
89
Poly-Si：Ar/Cl2プラズマ表面反応機構(3)
p0
SiCｌ2(s)+ Si(s) + Cl+(g)  SiCｌ2(g)+ SiCｌ(s)
SiCｌ3(s)+
Si(s)
Cl+(g) 
SiCｌ4(g)
+Cl2+(g)  SiCｌ2(g)
SiCｌ(s) + Si(s) + Cl2+(g)  SiCｌ2(g)+ SiCｌ(s)
SiCｌ2(s)+ Si(s) + Cl2+(g)  SiCｌ2(g)+ SiCｌ2 (s)
SiCｌ3(s)+ Si(s) + Cl2+(g)  SiCｌ4(g)+ SiCl (s)
0.16 ( E/E0 -1 )
0.19 ( E/E0 -1 )
1/2
1/2
0.13 ( E/E0 -1 )
1/2
E/E
0.16 (
0 -1 )
1/2
0.16 ( E/E0 -1 )
1/2
0.19 ( E/E0 -1 )
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
90
Si 形状 (Cl,Cl+,Ar+,Cl2+：1.6x106,5x104,1x104,1x104)
Si
T=200[s]
T=600[s]
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
T=1000[s]
91
入射粒子、反応生成粒子の密度
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
92
ボッシュプロセスによる
Si エッチング
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
93
ボッシュプロセス
SF6/C4F8 9[sec]/9[sec] 24回
流入サンプル数：エッチング ステップ
F, SiF_x(+): 3x106, 1.5x105
デポジション ステップ CxFy, I(+) : 1.5x106, 3x105
# エッチングステップ 反応
Si + 4F
→ SiF_4
Si + SiF_x(+) → SiF_x
Polymer + SiF_x(+) → SiF_x
Polymer + F
→ CxFy
# デポジションステップ 反応
Mask + CxFy → Mask + Polymer
Polymer + CxFy
→ Polymer + Polymer
Polymer + I(+)
→ CxFy + I
Si + CxFy → Si + Polymer
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
94
ボッシュプロセス
初期形状
ポリマー膜
Si 形状
Si
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
95
参考文献（追記）
1）菅井秀郎
プラズマエレクトロニクス
オーム社、2000
2）M. A. Lieberman and A. J. Lichtenberg
Principles of Plasma Discharges and Materials Processing, Second Edition
Wiley-Interscience, John and Wiley and Sons, Inc. 2005
3）Kenichi Nanbu and Masakazu Shiozawa
kinetic Modeling of Rarefied Plasmas and Gases in Materials Processing,
AIP Conference Proceedings Vol.663,
Rarefied Gas Dynamics: 23ｒｄ International Symposium、2003
4）ミシガン大：Kushner教授のHP
ハイブリッドモデル（静磁場有無）
http://uigelz.eecs.umich.edu/
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
96
おわりに
プラズマシミュレーションの課題
・電子衝突断面積、重粒子間反応レート、ラジカル種
の壁面での失活レートなどのデータ整備
・実測データによる検証
・表面波プラズマ
・計算時間の短縮
・３次元化
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
97
参考文献
[1]C.K. Birdsall and A.B. Langdon
Plasma Physics via Computer Simulation
McGraw-Hill, New York, 1985
[2] Kenichi Nanbu
Probability Theory of Electron-Molecule,Ion-Molecule,Molecule-Molecule
, and Coulomb Collisions for Particle Modeling of Materials Processing
Plasmas and Gases,
IEEE Trans. Plasma Sci. Vol.28, 971 (2000).
[3]G.A. Bird.
Molecular Gas Dynamics and the Direct Simulation of Gas Flow
Clarendon Press Oxford, 1994.
[4]南部健一
原子・分子モデルを用いる数値シミュレーション 第3章
モンテカルロ法の基礎. コロナ社, 1996.
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
98
[5]E.W. McDaniel and E.A. Mason:
The Mobility and Diffusion of Ions in Gases Wiley, New York, 1973
[6]K.Okazaki,T.Makabe and Y.Yamaguchi
Appl. Phys. Lett. Vol.54, 1742 (1989)
[7]T.Makabe, N.Nakano and Y.Yamaguchi
Phys. Rev. A Vol.45, 2520 (1992)
[8]T.J. Sommerer and M.J. Kushner.
Journal of Applied Physics, Vol.71, pp. 1654--1673, 1991.
[9]真壁利明
プラズマエレクトロニクス 培風館, 2000.
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
99
[10]H.M. Wu, D.B. Graves and M. Kilgores.
Plasma Source Sci. Technol., Vol.6, pp.231--239,1997.
[11]P.L.G. Ventzek, M. Grapperhaus, and M.J. Kushner.
J. Vac. Sci. Technol. B 12(6), pp. 3118--3137, 1994.
[12]M. Shiozawa and K. Nanbu.
Coupling of plasma and flow in materials processing.
Thin Solid Films Vol.457,pp.48--54,2004.
Copyright 2002-2008 PEGASUS Software Inc., All rights reserved.
100