第18回プラズマエレクトロニクス講習会プラズマプロセスモデリングと将来技術の最前線実用シミュレーション技術

第18回プラズマエレクトロニクス講習会プラズマプロセスモデリングと将来技術の最前線実用シミュレーション技術
第18回プラズマエレクトロニクス講習会プラズマプロセスモデリングと将来技術の最前線実用シミュレーション技術 2007年11月2日ペガサスソフトウェア株式会社松永史彦、中舘博 Copyright PEGASUS Software Inc., All rights reserved.

シミュレーションの必要性プロセス開発の現状
ブラックボックス装置の運転条件幾何形状、ガス種、ガス圧、パワー等　　　装置内プラズマやラジカルの状態　　　基板薄膜や微細構造の状態最終的な製品の状態を評価して、運転条件にフィードバックする。装置内の状態を測定するのは困難で、ブラックボックスとして扱われる。 Copyright PEGASUS Software Inc., All rights reserved.

シミュレーションの必要性シミュレーションの有効利用
基板に到達するイオン、ラジカル、スパッタ粒子のフラックスエネルギー分布入射角度分布　など数値シミュレーションにより装置内の状態を評価する。装置の運転条件幾何形状、ガス種、ガス圧、パワー等　　　装置内プラズマやラジカルの状態　　　基板薄膜や微細構造の状態装置内の状態を可視化することで、プロセス開発の効率を向上させることを目指す。 Copyright PEGASUS Software Inc., All rights reserved.

プラズマ装置内に生じている現象方程式の非線形性
ポアソン方程式荷電粒子の運動方程式荷電粒子の生成 Ra = （電子エネルギー分布の関数）荷電粒子の運動は電位分布（電界）の影響を受ける荷電粒子が運動すれば空間電荷分布が変化し、電位分布も変化する荷電粒子の生成率は電子エネルギー分布の関数であり、電子エネルギー分布は電位分布の関数である非線形性が非常に強い。短いタイムステップで時間を追跡していくしかない。 Copyright PEGASUS Software Inc., All rights reserved.

プラズマ装置内に生じている現象気相反応（CF4の場合）
Copyright PEGASUS Software Inc., All rights reserved.

Particle-In-Cell / Monte Carlo Collision Method プラズマの粒子シミュレーション法
PIC/MCC法 Particle-In-Cell / Monte Carlo Collision Method プラズマの粒子シミュレーション法 Copyright PEGASUS Software Inc., All rights reserved.

PIC/MCC法荷電粒子の壁面での境界条件
+ イオンが壁面に入射すると、中性化して反射する固体の内部に進入するイオンは壁面で消滅するとして取り扱う壁面　 Copyright PEGASUS Software Inc., All rights reserved.

イオンが壁面に入射すると、電子が放出される場合がある。入射者イオン1個あたり平均何個の電子が放出されるか→2次電子放出係数金属材料では　0.1前後の値をとるものが多い。MgO のように　3.0　という大きい値をとる材料もある。 DCプラズマの場合には2次電子放出がプラズマを維持するため、2次電子放出係数は特に重要。 + - 壁面　 Copyright PEGASUS Software Inc., All rights reserved.

希薄気体流れ場の粒子シミュレーション方法
DSMC法とテスト粒子モンテカルロ法希薄気体流れ場の粒子シミュレーション方法 Copyright PEGASUS Software Inc., All rights reserved.

DSMC法 Direct Simulation Monte Carlo Method 直接シミュレーションモンテカルロ法、モンテカルロ直接法 1960年代に G.A. Bird が最初に始めた。希薄気体の流れ解析に適する方法であるが、原理的には（計算量の問題を無視すれば）連続流の領域も解析可能。 Copyright PEGASUS Software Inc., All rights reserved.

DSMC法最大衝突数法(1) 粒子 i と粒子 j が衝突する確率は、i と j の相対速度の関数（剛体球の場合は比例）である。セル内の全ての2粒子の組み合わせの相対速度を評価するのは、計算コストがもったいない。例えば、セル内に粒子が 20個あれば、20x19=380対の相対速度を評価しなくてはならない。 Copyright PEGASUS Software Inc., All rights reserved.

DSMC法流入/流出境界(5) 粒子は流入境界と流出境界のいずれからも入ってくるし、また、出て行く。境界から出て行くフラックスは、計算結果として求まる量である。境界から入ってくるフラックスは境界条件として与える量である。 DSMC法では流入境界と流出境界に区別はない。粒子の出入りの正味のバランスの計算結果として、流入境界か流出境界かが決まる。 Copyright PEGASUS Software Inc., All rights reserved.

連続体流体と希薄気体の振舞の違いの例(1)
速度Vで壁が移動クエット流れ流体速度分布は？静止 Copyright PEGASUS Software Inc., All rights reserved.

2体衝突近似に基づくモンテカルロシミュレーション
Monte Carlo simulations using the binary collision approximation (MC-BCA) Copyright PEGASUS Software Inc., All rights reserved.

BC-MCA 計算例 Depth Profile

マグネトロンスパッタ装置２．PIC/MCC法プラズマ解析

マグネトロンスパッタ装置３．MC-BCA スパッタリング解析
スパッタリングフラックス分布（ターゲットの損耗分布）スパッタリング粒子放出角度分布 Copyright PEGASUS Software Inc., All rights reserved.

マグネトロンスパッタ装置４．スパッタ粒子の輸送

流体モデル前提平均自由行程　 λ ∝　1/p （p : ガス圧）装置の代表長さ　L P　大　　　λ　小　　　粒子モデルでは計算がきつくなる。 λ/L ≪ 1 なら電子・イオンの集団をそれぞれ連続体として取り扱うことができる。例　アルゴンプラズマの場合電子温度 1 [eV] として　λe～8.4×10-5/pTorr λe/L ≪1 pTorr≫8.4×10-5/L Ar+ 温度　室温として　 λi～4.8×10-5/pTorr λi/L ≪1 pTorr≫4.8×10-5/L p > 10 ～100 mTorr Copyright PEGASUS Software Inc., All rights reserved.

流体モデルマクロ量数密度： nα(r,t)=∫fα (r,cα,t)d3cα 流速： vα(r,t)=∫ cα fα (r,cα,t)d3cα/ nα(r,t)=< cα> エネルギー：<εα>(r,t)=< mαc2α/2>= ∫ (mαc2α/2) fα (r,cα,t)d3cα/ nα(r,t) r：位置座標, t：時刻, c：粒子の速度ベクトル m：粒子の質量 f(r,c,t)：速度分布関数添字α：粒子の種類(e,ion) Boltzmann 方程式これらの量をもちいて運動を記述する。 Copyright PEGASUS Software Inc., All rights reserved.

流体モデル基礎方程式のセット連続の式 ∂ne ∂t + ▽・Γe = Re ∂nion ∂t + ▽・Γion = Rion 運動方程式 Γe = - neμe E - De▽ne Γion = ± nionμion E - Dion▽nion 電子エネルギーバランス式 +eΓe・E = Reng (ne<εe>) +▽・ 5 3 Γe<εe> ∂ ∂t neDe▽<εe> - 電界 - ▽・( εrε0 ▽φ) = e ( ΣionnionZion – ne ) , E = - ▽φ Γα:粒子フラックス (≡nαvα) μα:移動度 Dα:拡散係数 Rα:粒子の生成レート E:電界 Reng:電子エネルギーの損失レート Copyright PEGASUS Software Inc., All rights reserved.

流体モデル計算手順これらの基礎方程式を離散化（差分法）して解くコイル電流から生ずる誘導電界・パワー電位に関するポアソン方程式 E これが一番時間がかかる Γe、Γion を計算モンテカルロ法による電子エネルギー分布関数の計算　ハイブリッドモ　デルの場合連続の式 ne, nion 電子エネルギーバランス式 <εe> この間、中性粒子種の密度は変化しないと考える中性粒子種に関する流れの計算 t=t+Δt 数周期毎に１回やり取りする Copyright PEGASUS Software Inc., All rights reserved.

流体モデル基礎方程式連続の式 ∂n ∂t + ▽・nv = R 運動方程式 + ▽・mnvv + ▽・P - neZE= R mom ∂mnv ∂t エネルギーバランスの式 ∂n<ε> + ▽・nv<ε> + ▽・P ・v + ▽・Q + env・E= R eng ∂t R:電離等の反応による粒子の生成レート Rmom:衝突による運動量の交換レート Reng:衝突によるエネルギーの交換レート P:応力テンソル=mn<(c-v)(c-v)>=mn<c’c’> Q:熱流ベクトル=mn<c’2c’>/2 Copyright PEGASUS Software Inc., All rights reserved.

流体モデル生成レート（電子の場合） νi　：電離の衝突周波数 νa　：電子付着の衝突周波数 νm ：運動量移行衝突周波数 νex　：励起の衝突周波数 Ei　：電離エネルギー Eex　：励起エネルギー N　　：衝突相手分子の数密度 M　　：衝突相手分子の質量　 σk　　：k衝突の衝突断面積 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 （分子は電子に対して静止しているという近似）・各種衝突について、衝突断面積のデータが必要・電子の速度分布関数を知る必要 Copyright PEGASUS Software Inc., All rights reserved.

流体モデル運動方程式の簡略化（ドリフトー拡散モデル）・応力テンソルPは等方的  スカラー圧力 p p=mn<c’2>/3 = 2nεT/3 = nkBT T:電子温度 + ▽・mnvv + ▽nkBT + neE = - v ∂mnv ∂t 1 mnνm 左辺第１項のオーダー : 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 PEGASUS Software Inc., All rights reserved.

流体モデル運動方程式の簡略化（ドリフトー拡散モデル）・τm/τv ≪ 1 の場合は第１項（時間微分）は無視できる。アルゴン中の電子の場合　νm ～ 1010ptorr [/s] τv=107 [s]（10MHzに相当）とすると τm/τv～10-3/ptorr 　　　　p > 10 [mTorr] なら条件が成り立つ。・v/vT ≪ 1 の場合は第２項（対流項）は無視できる。プロセスプラズマではだいたいこの条件は満たされていると考えてよい。 Copyright PEGASUS Software Inc., All rights reserved.

流体モデル運動方程式の簡略化（ドリフトー拡散モデル）圧力勾配項 ▽nekBTe meneνm 1 ne ▽( nevT2 3νm neτmvT2 3 ) = = 1 ne ▽(neDe)～ De ▽ne eE meνm = μeE 外力(電界）項 μe = e meνm 移動度 De = τmvT2 3 μekBTe e = kBTe meνm 拡散係数 Γe = - neμe E - De▽ne フラックス Copyright PEGASUS Software Inc., All rights reserved.

流体モデルエネルギーバランス（電子） ∂n<ε> + ▽・nv<ε> + ▽・P ・v + ▽・Q + env・E= R eng ∂t ・熱流ベクトルは温度勾配に比例 Q= - K▽Te = ▽εT = neDe▽εT 　　(K=5De/2) 2K 3kB 5 3 ・圧力はスカラー ▽・P・ve = ▽・(pve) = ▽・(neεTve) = ▽・(ΓeεT) 2 3 ・平均速度の運動エネルギーは、熱運動のエネルギーに対して無視 <ε> ≒ εT ≫ meve2/2 ・エネルギー損失 2me M Ei <ε> Eex Reng = -ne(νm + νi νex )<ε> = -neνε<ε> Copyright PEGASUS Software Inc., All rights reserved.

流体モデル電子に関する基礎方程式のセット ∂ne ∂t + ▽・Γe = (νi -νa) ne 連続の式運動方程式 Γe = - neμe E - De▽ne エネルギーバランス式 (ne<ε>) +▽・ 5 3 Γe<ε> ∂ ∂t neDe▽<ε> - +eΓe・E = -neνε<ε> ・νε≪ νm であるのでエネルギーバランス式においては　一般に時間微分項を無視することはできない。 Copyright PEGASUS Software Inc., All rights reserved.

流体モデルイオンに関する基礎方程式のセット・運動方程式に関しては、電子と同様ドリフトー拡散モデル・イオンのエネルギーは電離や電子付着にほとんど影響しない　ので、エネルギーバランス式は解かない。イオン温度は一定　と仮定。連続の式 ∂nion ∂t + ▽・Γion = νine ( for positive ion) νane ( for negative ion) 運動方程式 Γion = ± nionμion E - Dion▽nion μion = (mion+M)e|Z| mionMνm,ion Dion = μionkBTion e Copyright PEGASUS Software Inc., All rights reserved.

流体モデル電界についての基礎方程式・電子／イオンの運動を計算するために電界Eを求める必要・電気的に中性とはいえ、シースなどでは電気的中性　が破られる。( nion > ne ) ・電子／イオン密度の空間分布の変化により電界も変化電位に関するポアソン方程式を解く - ▽・( εrε0 ▽φ) = e ΣαnαZα E = - ▽φ φ: 電位　　　　 ε0 : 真空の誘電率 εr : 比誘電率　（プラズマ領域では 1） Copyright PEGASUS Software Inc., All rights reserved.

電子の速度分布関数に対する仮定衝突周波数： νk/N = ∫σk(c)f(c)c d3c （１）局所マクスウェル分布を仮定 4π
流体モデル電子の速度分布関数に対する仮定衝突周波数：　νk/N = ∫σk(c)f(c)c d3c （１）局所マクスウェル分布を仮定 4π c2 exp me 2πkBTe 3 2 mec2 2kBTe - f(c)=fLM(c;Te)= νk/N はTeの関数色々なTeの値に対してνk /N を計算し、Te～ νk /Nの関係をテーブルとして用意しておく。連続の式運動方程式（粒子フラックス） μe, De R, Rmom,Reng 電子エネルギーバランス式 Te（セル毎） νk 電位のポアソン方程式 Copyright PEGASUS Software Inc., All rights reserved.

流体モデル電子の速度分布関数に対する仮定（２）一様電界中の電子スウォームの分布関数を用いる一様な電界Eの中を移動する電子群に対する２項近似のボルツマン方程式を解き分布関数を求める。 Nμe、NDe、νk/N および<ε>はそれぞれ換算電界E/Nの関数色々なE/Nの値に対して計算を行い、 Nμe、NDe、νk/N ～E/N またはNμe、NDe、νk/N ～ <ε>関係のテーブルを作成しておく。連続の式運動方程式（粒子フラックス） μe, De R, Rmom,Reng E/Nの関数電位のポアソン方程式 E/N（セル毎） (局所電界近似) または μe, De R, Rmom,Reng 電子エネルギーバランス式 <ε>（セル毎） <ε>の関数 Copyright PEGASUS Software Inc., All rights reserved.

流体モデルハイブリッドモデル（粒子モデル＋連続体モデル）（３）電子のモンテカルロシミュレーションを用いるプラズマの連続体モデルの計算でえられた電子密度、分子密度、および(一般には時間依存の）電界E(r,t)をもちいて、同じ計算体系（形状）で電子のモンテカルロシミュレーションをおこなう。電子の速度分布関数f(r,c)をえる。（セル毎）（　実際にはエネルギー分布関数f(r,ε)　） μe, De R, Rmom,Reng νk　および Te がセル毎に求まる。連続の式 ne(r)（周期平均）　　電界E(r,t) 運動方程式（粒子フラックス）電位のポアソン方程式 Copyright PEGASUS Software Inc., All rights reserved.

流体モデルイオンの輸送パラメータ電子と異なる点・質量　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 PEGASUS Software Inc., All rights reserved.

流体モデルイオンの輸送パラメータ（１）　実測データを使用する・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 PEGASUS Software Inc., All rights reserved.

流体モデルイオンの輸送パラメータ（２）　次の近似式からμを求める。 Nμ= 3e|Z| 8 1 Ω(1,1)(Teff) π(mion+M) 2mionMkBTeff 1/2 kBTeff=kBTn +(1/3)Mvd2 Ω(1,1)(Teff)= (kBTeff)-3∫ε2σm(ε)exp(-ε/kBTeff)dε 1 2 ∞ D = μkBTeff e Copyright PEGASUS Software Inc., All rights reserved.

流体モデル電磁場とのカップリング ICPの場合 ▽2Eθ – iωμ0σEθ= 0　を解く Eθ σ|Eθ|2(Power) Eθ：誘導電界 ω：コイル電流の角振動数 μ0:真空の透磁率 e2ne me(νm+iω) σ：プラズマの導電率 = 連続の式 ne(r)（周期平均）　　 νm 電子モンテカルロ計算（ハイブリッドの場合）運動方程式電子エネルギーバランス式電位のポアソン方程式 Copyright PEGASUS Software Inc., All rights reserved.

流体モデル中性ガス流れとのカップリング原子、励起種生成衝突による励起、解離原料ガスの消費ガスの流動中性粒子種の密度変化電離、電子付着などに影響・中性ガスの流れシミュレーション　　各種成分の質量保存式（対流－拡散形）　　混合ガスの運動方程式　　混合ガスのエネルギー保存式 nα vgas （α:中性粒子種） Tgas ・電子・イオン運動の時間スケールとガス流れの時間　スケールはオーダーが異なるので、分離しておこなう。 Copyright PEGASUS Software Inc., All rights reserved.

流体モデル境界条件・境界面における粒子フラックスを規定する。 ^ ^ ^ (ξnμE・n) n （ξE・n > 0 のとき） Γdrift = ξ=±1 ^ （ξE・n > 0 のとき） Γdiff = (-D ▽n・n ) n ^ （n境界=0） Γout = aΓdrift + bΓdiff (a,b=0 または1) あるいは Γout = (1/4)nvT n , vT=√(8kBT/πm) ^ n: 境界面における外向き単位法ベクトル Copyright PEGASUS Software Inc., All rights reserved.

流体モデル境界条件通常、電子に関しては ^ Γ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 PEGASUS Software Inc., All rights reserved.

流体モデル離散化 nij μij Dij Teij φij Γr i+1/2,j Γr i-1/2,j Γz i,j+1/2 Γz i,j-1/2 Er i+1/2,j Er i-1/2,j Ez i,j-1/2 Ez i,j+1/2 セル(i,j) ・２次元直交メッシュ軸対称平面２次元・変数配置　　　staggered ・コントロールボリューム法 Copyright PEGASUS Software Inc., All rights reserved.

流体モデル Sharfetter – Gummel スキーム(指数法) vi+1/2 ni ni+1 Γi+1/2 Δx Γi+1/2 = vi+1/2 exp(αΔx)ni - ni+1 exp(αΔx) - 1 α= (v/D)i+1/2 = (±μED)i+1/2 ・Γ=nv-D(dn/dx) として、Γi+1/2を評価するのに、１次元定常問題　dΓ/dx=0 の解析解を局所的に使用。・|α|Δx ≫1 では１次の風上法、 |α|Δx ≪1では中心差分になる。計算が安定 Copyright PEGASUS Software Inc., All rights reserved.

流体モデルポアソン方程式のsemi-implicit 解法 ∂ρ ∂t = - ▽・Σe ZαΓα アルゴン　 p=100 [mTorr] ne=1012 [/cm3] = - ▽・Σe |Zα|nαμαE ～ - e(Σ|Zα| nαμα) ▽・ E ρ ε0 τM=1.5×10-12 [s] = - e(Σ|Zα | nαμα) eΣ|Zα| nαμα ε0 誘電緩和時間:　τM = 陽解法をもちいる場合　Δt < τM とする必要 Copyright PEGASUS Software Inc., All rights reserved.

流体モデルポアソン方程式のsemi-implicit 解法 - ▽・(ε0 ▽φ) = e ΣαnαZα nα= nα,old + Δt (-▽・Γα+Rα) = nα,old - Δt (|Zα|/Zα)μα nα ▽・E + Δt (Dα ▽2nα,old +Rα) 1 ε0 eΣαZα nα,old + Δt (-▽・Γα,diff + Rα) Δt τM 1+ ▽2φ= - Δt > τM としても安定に計算できる。 Copyright PEGASUS Software Inc., All rights reserved.

流体モデル時間刻みΔt Δt < min Δr/(μeEr), Δz/(μeEz) ・電界ドリフトによる移動が１ステップで１メッシュを　こえない。・電子のドリフト速度で評価・生成レートが大きい場合には、この評価よりもかなり　小さくしなければならない場合も有る。・だいたいΔt =10-10～10-12 程度 Copyright PEGASUS Software Inc., All rights reserved.

ドライエッチング装置反応式 e+F → F(+)+2e e+F2 → e+F2 e+F2 → F(-)+F e+F2 → F+F+e e+SF6 → e+SF6 e+SF6 → SF6 (-) e+SF6 → SF5(+)+F+2e e+SF6 → SF4(+)+2F+2e e+SF6 → SF3(+)+3F+2e e+SF6 → SF2(+)+4F+2e e+SF6 → SF(+)+5F+2e e+SF6 → SF5(-)+F e+SF6 → SF3+3F+e e+SF6 → SF2+4F+e e+SF6 → SF5+F(-) e+SF5 → SF4(+)+F+2e e+SF5 → SF5(+)+2e e+SF4 → SF4(+)+2e e+SF3 → SF3(+)+2e e+SF2 → SF2(+)+2e e+SF → SF (+)+2e X(+)+Y(-) →X+Y (X,Y∈SFn,F) 壁面反応すべての壁面 F → 1/2F2 F2, SFn(n=1,5)は壁面で0.1%損失 Copyright PEGASUS Software Inc., All rights reserved.

おわりに・熱、構造、流体、電磁気分野ではシミュレーションが定着・プラズマ分野ではシミュレーションの普及はまだこれから　　　　　・計算時間の短縮　　　　　・断面積等のデータ整備　　　　　・実験とシミュレーション結果の照らし合わせ　　　　　・３次元化　　　　　・使い勝手　　　　　・メーカー・大学等研究機関・フトウェア開発者　　　　　　の連携装置設計により役立つシミュレーションツールを目指す Copyright PEGASUS Software Inc., All rights reserved.

第18回プラズマエレクトロニクス講習会プラズマプロセスモデリングと将来技術の最前線実用シミュレーション技術

Similar presentations

Presentation on theme: "第18回プラズマエレクトロニクス講習会プラズマプロセスモデリングと将来技術の最前線実用シミュレーション技術"— Presentation transcript:

Similar presentations

About project

フィードバック

ログインする

Auth with social network:

第18回プラズマエレクトロニクス講習会 プラズマプロセスモデリングと将来技術の最前線 実用シミュレーション技術

Similar presentations

Presentation on theme: "第18回プラズマエレクトロニクス講習会 プラズマプロセスモデリングと将来技術の最前線 実用シミュレーション技術"— Presentation transcript:

Similar presentations

About project

フィードバック

第18回プラズマエレクトロニクス講習会プラズマプロセスモデリングと将来技術の最前線実用シミュレーション技術

Presentation on theme: "第18回プラズマエレクトロニクス講習会プラズマプロセスモデリングと将来技術の最前線実用シミュレーション技術"— Presentation transcript: