中高圧力下のスパッタリングにおける粒子輸送過程
Particle transport process in high pressure sputter deposition
中野 武雄∗

Takeo NAKANO
成蹊大学工学部
Dept. of Applied Physics, SEIKEI Univ.

1

はじめに低圧放電プラズマを用いたスパッタ製膜プロセスでは、基板面での膜厚分布や飛来する粒子のエネルギー分布・角度分布を理論的に予測することが難しい。ひとつには、ターゲットから飛び出した粒子が、放電用に導入されている希ガスによって散乱されて方向を変えるからである。このため、真空蒸着で用いられているような幾何学的な計算が使えない。またもうひとつには、ターゲットから飛び出した直後の粒子の運動エネルギー (数 eV程度) が、雰囲気の希ガスの熱運動エネルギーに比べて非常に大きいからである。このため、粒子の輸送過程を単なる拡散とみなすこともできない。
この「ターゲットから放出された高エネルギーの粒子の輸送過程」を扱う手法として、
モンテカルロ法 (MC 法) を用いた計算機シミュレーションが 1980 年代の前半から提案されてきた。これは、主にガスとの衝突からなるスパッタ粒子の輸送過程を Markov 連鎖で取り扱い、10 万 ∼ 100 万程度の粒子を追跡することで、壁面へ到達する粒子の位置分布やエネルギー・角度の分布を得ようとするものである。1983 年からは元廣らがスパッタ粒子の減速過程を扱うために [1]、また 1984 年には Somekh がターゲット起源の高エネルギー粒子によって成長中の薄膜が受けるダメージに関する興味から [2]、それぞれ MC法によるシミュレーションモデルを提案した。1989 年にはスパッタ粒子の平均自由行程がその速度に依存することを Turner らが指摘し、さらに衝突するガスにも Maxwell 分布の熱運動を仮定して、MC シミュレーションの精度を向上させた [3]。このようにモデルが確立されたため、1990 年代に入ってからは、MC 法を利用したスパッタ製膜プロセスの研究が数多く報告されるようになった。昨今では、計算機の急激な高性能化に伴い、このようなシミュレーションを手元のワークステーションや PC で行うことも充分可能になってきている。これらの研究例についてや、シミュレーションモデルの解説としては、
元廣による総説記事 [4] が参考になる。
我々は、もともとは LaB6 スパッタ膜の組成が圧力に依存して変化する現象を説明するために、MC 法による計算を扱いはじめた。その後、ガスの熱運動をより厳密に取り扱う方法 [5] や、衝突によって速度が落ちてきたスパッタ粒子の輸送過程を拡散的に取り扱∗〒 180-8633 東京都武蔵野市吉祥寺北町 3-3-1, E-mail: nakano@apm.seikei.ac.jp

い、これを MC 法と組み合わせる方法 [6] を提案してきた。これらによって、原理的にはスパッタ粒子の輸送過程を 10 Pa 以上の高い圧力でも取り扱うことが可能になった。
ガスの圧力が高くなると、いわゆる「熱中性化」するスパッタ粒子の割合が多くなる。
高圧力下では平均自由行程が短くなるので、スパッタ粒子は気相中で多数回衝突し、当初の大きな運動エネルギーを失っていく。このようにして、スパッタ粒子の運動エネルギーが雰囲気ガスの熱エネルギーと同程度にまで減速されることを熱中性化 (thermalization)と呼ぶ。熱中性化の前には、スパッタ粒子はターゲットから離れる向きの運動量を持っており、あたかも静止したガスにあたって飛び散る弾丸のように (ballistic に) 輸送される。
逆に熱中性化後のスパッタ粒子には、ガス原子が四方八方から当たることになり、スパッタ粒子は拡散的 (diffusive) に輸送される。中〜高圧力でのスパッタ製膜においては、このふたつの過程の競合が、スパッタ膜の膜厚分布や、化合物スパッタ膜の膜組成などを支配する。
本稿では、この熱中性化前後の輸送過程の競合過程を、MC シミュレーションによって再現した例を紹介する。また我々の導入した手法を用いると、熱中性化後も各々のスパッタ粒子の運動を追跡でき、スパッタ粒子の気相中の空間密度が評価できることも示す。

2

シミュレーションの方法

図 1 に、我々が用いている MC シミュレーションの流れ図を示す。シミュレーション
中のスパッタ粒子は、角度とエネルギー (つまりベクトル的な意味での「速度」)、および位置を内部変数として持っている。
まずスパッタ粒子はターゲット上から、ある角度・エネルギーで放出される。続いてある自由行程だけ飛行し、ガス原子と衝突する。この途中の経路に容器の壁面があると、スパッタ粒子はそこに付着する。ガスとの衝突の際には、ガスの速度も計算され、二体間のポテンシャルによって散乱プロセスが処理される。新しい角度とエネルギーとが求まると、再び次の自由行程が繰り返される。衝突後のスパッタ粒子のエネルギーがあらかじめ決めておいた閾値 (ガスの熱エネルギー等) を下まわった場合には、スパッタ粒子はそこで熱中性化したとみなされ、その位置が記録されて追跡は打ち切られる。
試行数分の追跡が終わると、熱中性化したスパッタ粒子の分布から、これらが壁面のどこに拡散付着するかが計算される。この拡散付着した粒子の分布と、MC 法での追跡中に付着した粒子の分布とが加えられ、最終的な分布を得る。
以降の節では、それぞれのプロセスをもう少し詳しく説明する。

2.1

ターゲットからの放出

スパッタ粒子のターゲットからの放出位置は、実際に製膜実験を行った後のターゲットに生じるエロージョントラックの深さ分布から、損失体積に比例するように決定する。例えば円形ターゲットの動径方向の深さ分布を d(r) とすれば、放出位置の動径成分は 2πrd(r)に比例する。
放出の際のエネルギー分布としては、Thompson による式 [7] を用いるのが一般的である。これはエネルギー E1 を持った一次イオンによってターゲット内部に引き起こされたカスケードが、エネルギーを減衰させながら拡散していく過程を扱ったもので、粒子の放図 1: モンテカルロシミュレーションの流れ
出エネルギー E の分布関数 f (E) は

f (E) ∝ E

1−q(Eb + E)/ΛE1(E + Eb )3
Λ = 4mi mt /(mi + mt )2(1)で与えられる。ここで Eb はターゲット材料の原子結合エネルギー、mi と mt はそれぞれイオンとターゲット材料原子の質量である。
Thompson の理論によれば、放出角度分布は余弦則に従うことになる。すなわち放出角度の法線からの角度を θ とすれば、分布関数は cos θ に比例する。放出角度の分布については、余弦則に比べて斜め方向の分布が大きな蝶ねじ型の分布になる、という報告も多い [8] が、我々のシミュレーションでは取り扱いが単純な余弦則を仮定している。

2.2

自由行程

いまガス密度を ng 、衝突断面積を σ とする。スパッタ粒子が充分速く、ガスの運動が無視できる場合には、スパッタ粒子の平均自由行程 λ は 1/(ng σ) で与えられる。逆にスパッタ粒子がガスよりも遅い極限を考えると、ガス原子が周囲からぶつかってくることになるので、衝突周波数 ν の方が一定に近づき、λ(= vp /ν) はスパッタ粒子の速度 vp に比例する。質量 mg のガス原子が温度 Tg で Maxwell 分布をしている場合の式は、前述のTurner の論文で次のように与えられている。
√πmg vp 21µ q¶
λ =mg2σng kTg Ψ vp 2kTg
Ψ(x) = x exp (−x2 ) + (2x2 + 1)

(2)Z x0exp (−y 2 )dy

平均自由行程 λ が与えられたとき、各々のスパッタ粒子の自由行程 l は、λ を平均値とするような Poisson 分布に従う。すなわち l の分布関数 f (l) ∝ exp(−l/λ) となる。

2.3

衝突・散乱過程

スパッタ粒子に衝突してくるガスの速度分布は、Maxwell 分布によって与えられる速度空間上での密度分布に、スパッタ粒子とそのガスとの相対速度 vR の重みを乗じたものになる。いまスパッタ粒子が z 軸に沿って飛行しているとして、球面座標系でとった速度空間中の微小要素 (vg , θ, φ) 〜(vg + dvg , θ + dθ, φ + dφ) からの衝突周波数への寄与 dν はsin θdθdφ
4π22= vp + vg − 2vp vg cos θdν = σvR ng fM (vg )dvgvR 2(3)で与えられる。fM は Maxwell の速さ分布である。この dν を全速度空間で積分したのが、前出の式 2 になる。
式 3 は衝突してくるガスの速度分布とみなせるので、これを用いれば MC 法で確率的にガスの速度を決定できる。この分布は vR の項を含むので、スパッタ粒子の速度 vp に依存する。
スパッタ粒子とガスの速度とがそれぞれ決まったら、二体間のポテンシャルによる散乱を計算する。ポテンシャル U (r) としては、中心対称の Born-Mayer 型斥力ポテンシャル

U (R) = A exp (−br)(4)

を採用した。ここで A と b は定数で、各原子のパラメータは Abrahamson によって与えられている [9]。弾性散乱を仮定すれば、スパッタ粒子–ガス間の相対速度 vR の絶対値は変化せず、向きだけが変わる。この偏向角度 θR は、重心系において

θR = π − 2Z ∞bc /r2qr0

1 − (bc /r)2 − 2U (r)/µv0 2dr(5)µ = mp mg /(mp + mg )で与えられる。bc は衝突パラメータ、v0 は無限遠における vR の値、r0 は角運動量保存則から求まる二体間の再接近距離である。mp と mg はスパッタ粒子とガスの質量である。
このようにして得られた衝突後の相対速度から衝突後のスパッタ粒子の速度を求め、座標系を重心系から実験室系に戻してやれば、衝突後のスパッタ粒子の方向とエネルギーとが決定できる。

2.4

拡散過程

圧力が増加すると衝突イベントが増えるため、MC 法で粒子が壁面に到着するまでを追跡すると、圧力の 2 乗に比例した計算時間がかかってしまう。そこで熱中性化以降の輸送は拡散方程式で取り扱った。
スパッタ粒子が熱中性化した後の拡散は、定常状態の拡散方程式

0 = D∆np (x) + ρ(x)(6)によって記述できるものと考える。ここで D は拡散係数、np はスパッタ粒子の空間密度、x は空間座標である。このとき壁面境界における粒子のフラックス F は、D grad npの法線成分によって与えられる。
ρ は粒子の湧き出しに対応する項であり、モンテカルロ法を用いた追跡中に熱中性化した粒子の座標 xTi h を用いてρ(x) =´
1 X 3 ³ Thδ xi − xTM C(7)iと表せる。ここで TM C は、モンテカルロ計算の試行数分のスパッタ粒子放出に対応する、実際のプロセスにおける時間である。このとき TM C 間に壁面に付着する粒子の面密度はTM C F となる。
いま η ≡ TM C Dnp なる量を定義すれば、式 6 は

0 = ∆η +X³

δ 3 xTi h´(8)iとなり、TM C 間の壁面への付着数面密度は grad η の法線成分で与えられる。つまり、拡散係数の値がわからなくても、拡散して壁面へ到達する粒子の分布を求めることができる。

3シミュレーションの適用例

3.1

基板ホルダ裏面へのまわり込み付着

スパッタ製膜では、スパッタ粒子のガスによる散乱のため、ターゲットから直接見えない容器壁面へも原子がまわり込んで膜が堆積する。元廣らは、このような背面へのまわり込みを、裏面/表面の膜厚比 (Back-Front ratio, BFR) で整理している [10]。彼らは BFRの圧力依存性を実験で求め、ある圧力で BFR が極大を持ち、高い圧力では減少することを示している。ただしこの頃には高圧力での信頼性の高いシミュレーションが (計算機能力などによる制限のため) 行えなかったため、高い圧力での減少については十分に議論されていなかった。
図 2 に、我々の装置で作成した Cu 膜での BFR の圧力依存性を示した。黒塗りの印が実験結果、中抜きの印がシミュレーションの結果である。製膜は図 2 の左に示したような円筒対称の装置を用い、DC マグネトロンスパッタで行った。半径 40 mm の Cu ターゲットが下部に置かれており、エロージョントラックの位置はほぼ半径 20 mm を中心に分布している。上部にある基板ホルダの半径は 40 mm で、その端から 5 mm 内側 (すなわち r = 35 mm) と 10 mm 内側 (r = 30 mm) とで BFR を求めた。r = 35 mm の方が図 2: 銅スパッタ膜の背面へのまわり込みホルダの端に近いため、BFR は大きな値になっている。背面での膜厚が薄いため実験値のエラーバーが大きく、シミュレーションと正確に一致しているかどうかを議論するにはやや不満が残るが、おおよその値と、高圧になると減少して一定値に近づく傾向は一致している。シミュレーションの結果に見られる、低圧から上昇してきてピークをとり、高圧に向けて緩やかに減少する傾向は、元廣らの実験結果とも一致している。
このような BFR の振る舞いは、スパッタ原子の輸送過程が ballistic なものから diffusiveなものに移行することが原因である。低圧の極限では、スパッタ粒子は直進するだけで、背面へは到達できない。ここから圧力が上昇すると、スパッタ原子はガスによる散乱を受けて方向が変わり、基板ホルダの背面に到達するようになってくる。この割合は最初圧力上昇とともに増加する。さらに圧力が増加すると、スパッタされた原子はある程度ターゲットから離れた場所で熱中性化するようになる。これらの原子は、その地点から拡散過程 (実際には平均自由行程スケールでのランダムウォーク) によって移動し、やはりある確率でホルダー背面に付着する。この段階になると、圧力の増加につれて熱中性化する地点がターゲット近傍に近寄ってくる。これは拡散原子の湧き出しがターゲットに近づくことを意味する。したがって基板ホルダー背面までの距離が相対的に遠くなり、まわり込みは減少する。高圧の極限では、粒子は放出後すぐに熱中性化し、輸送は拡散過程によってのみ支配される。すなわちターゲットのエロージョントラックが拡散粒子の湧き出しそのものとみなせるため、BFR はこの拡散方程式の解で決まることになり、一定値に近づく。
図 3 は、シミュレーションにおける熱中性化地点の分布が圧力によってどのように変化するかを示したグラフである。熱中性化した粒子の空間密度を濃度マップで表しており、図の上部に濃度スケールをリニアスケールの任意目盛で示してある。これはシミュレーションで 10 万個の Cu 原子をターゲットから放出させた結果である。熱中性化のエネル
ギー閾値は kTg としてある。グラフを見てわかるとおり、シミュレーションで BFR が最大値をとった 1.5〜2 Pa において、粒子の熱中性化が顕著になり始める。その後は圧力の増加とともに熱中性化地点の分布がターゲットに近づく結果がよく分かる。このため基板背面へのまわり込みの割合が一定に近づくわけである。

図 3: 銅スパッタ粒子の熱中性化位置の分布3.2化合物スパッタ膜の組成変化

合金や化合物など、多元素の物質をスパッタ製膜すると、膜の組成がターゲットとずれる現象が良く観察される。このずれは一般に放電ガスの圧力に依存する [11]。スパッタ粒子がガスによって散乱される過程はスパッタ粒子の質量に依存するので、成分元素ごとに輸送過程は異なり、膜組成にずれを生じさせうる。ここでは成分元素の質量が大きく異なる LaB6 (B/La〜11/139) を対象に、この圧力依存性を議論する。
図 4 に LaB6 スパッタ膜の組成の圧力依存性を示した。(a) は実際に行った製膜実験の結果、(b) がシミュレーションの結果である。化合物ターゲットでは 2.1 節で示した
Thompson のモデルを使えないので、ここでは B と La の放出時のエネルギーが等しいと仮定し、5 eV と 10 eV の場合を示した。B/La の値は、中心からの距離が異なる 3 つのトラック (r = 0 ∼ 5 mm, 10 ∼ 15 mm, 20 ∼ 25 mm) 上で求めたが、位置による依存性はあまり見られなかった。低圧力から圧力が増加すると B/La 比がいったん減少し、あるところから再び増加に転じるという、実験で見られた U 型の圧力依存性がシミュレーションでも再現されている。B/La の極小値はいずれも 4 程度と良い一致を見せている。
ただし極小となる圧力は、実験値が 10 Pa 前後であるのに対し、シミュレーションでは2〜3 Pa あたりの、より低い圧力となった。
なぜ B/La 比が図 4 に見られるような U 型の圧力依存性を示すのかを、B および La各原子の気相中の輸送過程の側面から考えてみる。スパッタされた B, La の各原子は、当初は数 eV という高いエネルギーを持っている。このときアルゴンガス原子の動きを無視すると、ボロン原子は アルゴン原子によって後方散乱されうるのに対し、ランタン原子では前方散乱しか起こらない。圧力が低く平均自由行程が長いと、スパッタ原子はあまりガス原子と衝突しない。したがって粒子はスパッタ直後の高いエネルギーを保ったまま、ballistic な輸送プロセスに支配される。この領域で圧力が上昇すると、衝突によって方向が大きく変えられやすいボロンの方が、基板に到達する確率がより大きく減少するため、膜におけるボロン組成が減少する。
一方、ある圧力以上になると B/La 比が化学量論比である 6 に向かって回復するのは、

図 4: LaB6 スパッタ膜の組成の圧力変化
スパッタ原子が熱中性化するからである。前節でも述べたように、圧力が上がると壁面に到達する前に空間中で熱中性化するスパッタ原子の割合が増える。また 2.4 節で述べたように、熱中性化した後の拡散段階において、定常状態での壁面への原子フラックスに拡散係数は影響しない。したがって原子の質量の違いによらず、同じ取り扱いができる。ただし熱中性化地点は La の方が B よりターゲットから遠くに分布するので、基板での原子フラックスは La の方が大きくなる。しかし高圧力での極限では、熱中性化地点はどんどんターゲットに近づき、エロージョントラックが拡散原子の湧き出しそのものになっていくから、基板への到達割合も B と La とで変わらなくなっていくと考えられる。
B/La 比が極小値をとるガス圧力が実験とシミュレーションとでずれている理由は、粒子輸送の立場から考えると、ballistic な輸送から diffusive な輸送に移行する圧力が異なるせいである。今回のシミュレーションでは、ガス圧力はチャンバーの全領域で一定と仮定していたが、高圧力・高放電電力におけるスパッタ製膜においてはターゲット近傍でガス密度が希薄化する効果が提唱されており [12]、このため実験での極小がより高い圧力になったのかもしれない。またシミュレーションではターゲットからの放出時のエネルギーを大きくすると極小の圧力が高くなったが、これは熱中性化により多くの衝突が必要となったため、平均自由行程の短い高圧力側にシフトしたと考えられる。本節のシミュレー
ションでは B 原子と La 原子のターゲットからの飛び出しエネルギーが等しく、かつすべての原子がこの一定のエネルギーで飛び出すとしていたが、この点についてより具体的なモデルが開発できれば、より定量的な議論が可能になるだろう。

3.3

スパッタ粒子の空間密度

最近、スパッタプロセスにおける、粒子の気相中での振舞いが注目を集めつつある。これはプロセスそのものに対する興味もさることながら、気相中でスパッタ粒子をイオン化する IPVD [13] などの手法が良い膜を得るのに有効であることがわかってきたからである。

図 5: モンテカルロ法で計算したスパッタ粒子の空間密度
スパッタ粒子の空間密度は、ターゲットからの粒子の放出頻度と、各粒子の滞在時間とから導くことができる。シミュレーションでこの両者を求めることを考えてみる。放出頻度は 2.4 節で出てきた TM C がわかれば、試行粒子数をこれで割ることによって求められる。さらに製膜実験で得られた基板上での成膜速度と、シミュレーションで得られた同じ場所での付着粒子面密度とを比べれば、TM C を求めることができる。滞在時間は、粒子の位置と速度を (熱中性化を無視して) ひたすら壁面に付着するまで追跡すれば計算できる。
このようにして、実験で得られた成膜速度と、シミュレーションでの付着原子面密度・滞在時間とから求めた、Cu 原子の空間密度プロファイルを示したのが図 5 である。2〜20 Pa の、4 つの異なる圧力での結果を示した。黒い線は容器の境界で、左側の一点鎖線はターゲット中心を通る対称軸である。スパッタ原子の放出は下側のターゲット面上の半径 2 cm あたりを中心に生じている。密度は 1016.5 ∼ 1019 m−3 の範囲を対数のカラースケールで示した。
Ar ガスの圧力が圧力が増加するとともに、原子密度も増加している。高い圧力ではスパッタ原子の平均自由行程が短くなり、拡散係数が小さくなる。これが高い密度をもたらしているのであろう。それぞれの圧力において、原子密度は粒子の湧き出しであるターゲットのエロージョントラックの近辺で最大となっている。一方、その他の容器境界では、近づくに連れて原子密度は減少している。シミュレーションでは、境界に到達したスパッタ原子は確率 1 でその場に付着するとした。これは拡散方程式では、境界での密度が 0 であることに対応する。容器壁面近辺での密度減少が高圧力下でより顕著になっているのは、原子が ballistic ではなくより diffusive に輸送されているからであろう。例えば低圧力の極限を考えると、ターゲットから飛び出した原子は全くガス原子と衝突せずに容器壁面までまっすぐ飛んでいく。この場合の原子密度プロファイルは、ターゲットからの原子フラックスが広がることのみから導かれ、容器壁面付近での急激な減少は示さないはずである。
今回得られた原子密度 1017 ∼ 1018 m−3 という値は、Dony らが吸光実験で得た値 [14] とほぼオーダー的に等しい。また我々が以前に行った、この銅プラズマの発光分析結果 [15]

も、これらの原子密度から半定量的に説明することが可能となっている。

4

まとめ

スパッタ粒子の輸送過程をモンテカルロ法でシミュレートする手法を紹介し、これを製膜プロセスにおける膜厚分布や化合物膜における組成の問題に適用した。中〜高圧力における製膜においては、数 Pa 程度を境にして起こる、ballistic な輸送と diffusive な輸送との競合によって、裏面へのまわり込みが極大となったり、化合物での軽元素の組成が極小値となることを示した。また MC 法を援用して、スパッタ粒子の空間プロファイルを評価する方法も紹介した。
高速粒子のガス中での減速・拡散という点では、現在のモデル化はほぼ満足できるものと我々は考えている。今後実験との一致を向上させるためには、ガス自体の密度・温度分布や、ターゲットからスパッタ原子が放出される過程について、既存のモデルとうまく組み合わせて行くことが重要となるだろう。

参考文献[1] T. Motohiro and Y. Taga, Surf. Sci. 134, L494 (1983).
[2] R. E. Somekh, J. Vac. Sci. Technol. 2, 1285 (1984).
[3] G. M. Turner, I. S. Falconer, B. W. James, and D. R. McKenzie, J. Appl. Phys. 65,3671 (1989).
[4] 元廣 友美, セラミックス 33(10), 805 .
[5] T. Nakano, I. Mori, and S. Baba, Appl. Surf. Sci. 113/114, 642 (1997).
[6] T. Nakano and S. Baba, Thin Solid Filims 343-344, 24 (1999).
[7] M. W. Thompson, Philos. Mag. 18, 377 (1968).
[8] 金原 粲. スパタリング現象. 東京大学出版会, (1984). 2.5 節.
[9] A. A. Abrahamson, Phys. Rev. 178, 76 (1969).
[10] T. Motohiro, J. Vac. Sci. Technol. A 4(2), 189 (1986).
[11] S. M. Rossnagel, I. Yang, and J. J. Cuomo, Thin Solid Films 199, 59 (1991).
[12] S. M. Rossnagel, J. Vac. Sci. Technol. A 6(1), 19 (1988).
[13] S. M. Rossnagel, J. Vac. Sci. Technol. B 16, 2585 (1998).
[14] M. F. Dony, J. P. Dauchot, M. Wautelet, M. Hecq, and A. Richard, J. Vac. Sci.
Technol. A 18(3), 809 (2000).
[15] Takeo Nakano, Noriyuki Ohnuki, and Shigeru Baba, Vacuum 59, 581 (2000).