中高圧力下のスパッタリング
における粒子輸送過程
中野武雄
成蹊大学工学部
2001年12月12日

内容
序論・背景
モンテカルロシミュレーション


1.
2.
3.

ターゲットからの放出過程
雰囲気ガスとの衝突
拡散過程の取扱い

粒子輸送過程の圧力依存性


1.
2.
3.

銅スパッタ膜の膜厚分布化合物スパッタ膜の組成変化（スパッタ粒子の空間密度評価）

スパッタ製膜プロセス(1)

Ar

C

→PUMP

スパッタ製膜プロセス(2)装置
・容器形状
・ガス (元素・圧力)・電力・マグネトロン磁場

プラズマ
・イオン生成・電場・磁場

ターゲット・エロージョン
トラックの形成
・材料原子(分子)の飛出し (θ, φ, E)・入射イオンの反跳

粒子輸送
・雰囲気ガスとの衝突・散乱
・拡散・(イオン化)

基板・表面拡散・ミキシング
・欠陥導入
・再スパッタ壁面:
粒子 flux
入射角分布
エネルギー

膜の性質
膜厚分布

粒子輸送過程 スパッタ直後(EP =1~10 eV)
 Ballistic な過程

2.散乱(EP>> kTg)
 Diffusive な過程(EP ~ kTg)

3. 熱中性化・拡散
1. 飛出しモンテカルロ・シミュレーション 1980年代前半に、Somekh, Motohiro らによって独立に開発された。
 スパッタされた粒子の “life” を構成する各過程を確率的に取り扱う。
 ターゲットからの飛び出し ガスとの衝突・散乱
 熱中性化後の拡散

 数十万〜数百万個の粒子を追跡し、膜厚分布や基板面への入射エネルギー分布などを計算する。

[x10-3]

モンテカルロ法の原理
3.0f(v)

2.0
1.000
1.05001000500

1000

F(v)

0.8
0.6
0.4
0.200

ある変数 v が適当な規則に従ってランダムに選ばれる場合、その確率密度関数 f (v ) を物理的なモデルに則って定義する。
 f (v ) に対する母関数 F(v )を計算する。
 計算機で 0~1 の乱数を発生させ、母関数から変数 vを求める。
v

シミュレーションの流れ
START
処理の流れ
ターゲットから放出
位置・角度・Energy

データの流れ自由行程の計算
位置壁面付着?

粒子付着位置の分布

Y

熱中性化位置の分布

N

拡散過程

散乱過程
角度・Energy

熱中性化?
NEND

Y
粒子数分終了?
NYターゲットからの放出(1)エネルギー分布Thompsonの式

Λ = 4 M I M T /( M I + M T ) 2
MI , MT: 入射イオン・ターゲット
原子の質量
Eb:
ターゲット材料の結合
エネルギー
E1:
入射イオンのエネルギー

1.0

dN/dE (a.u.)

1 − ( Eb + E ) / ΛE1f (E) ∝ E( E + Eb ) 3

Ar→CuEb = 3.5 eV
E1 = 400 eV

0.50.002040

Ejection Energy (eV)60

Depth (a.u.)

ターゲットからの放出(2)放出位置・角度f (θ ) ∝ cos θ0.0-0.5-1.001234Radial Position (cm)
4Y (cm)2

θ
0-2-4-4-20X (cm)24雰囲気ガスとの衝突
 ガスの速度分布には一様・等方な

Maxwell 分布を仮定
 粒子との相対速度の重みをかけて積分し、衝突周波数（および平均自由行程）を計算
 被積分関数を確率分布関数と考え、衝突してくるガスの速度とガス−スパッタ粒子間の相対速度とを決定する

衝突周波数sin θ dθ dϕdν = σ vR ng f (v g ) dv g
4πσ : 衝突断面積vRvg φf (v g ) : Maxwell分布vR = v p + v g − 2v p v g cos θ222vp

θ衝突周波数 (2)∞ν= ∫π2πdν∫∫θ ϕ⇐ (v g , θ ) → (v g , v R )vg =0 =0 =0∞= 2 π σ ngα 3 / 2 ∫ dv g
02σ ng kTg=
Ψ (αv p )π mg v p

λ=vpνvg
vpv p +vg

exp(−αv g )2∫ vR dvR
2v p −vg⇐ Ψ ( x) = x exp(− x 2 ) + (2 x 2 + 1)

平均自由行程π2Erf( x)自由行程
Mean Free Path (m )10v p >> v g-2

平均自由行程が一定v p << v g10衝突周波数が一定→MFP は粒子速度に比例-3Ar 400K, 1Pab max = 0.6 nm
10 -4
10100

1000

10000

Particle Speed (m /s)

毎回の自由行程 l: l
f (l) = exp − λ λ
1

衝突ガス速度の決定∞

ν = 2 π σ ngα 3 / 2 ∫ dv g
0vgv p +vg

exp(−αv g )2∫ vR dvR
2vpv p −vg
1444424444
3ν ' (v p , v g )Vg : w1 = ν ' (v p , Vg ) ν ' (v p , ∞)VR : w2 =VR∫ vR dvR
2v p −V gv p +V g

∫ vR dvR
2v p −V g

w1, w2: [0, 1) 間の一様擬似乱数

衝突ガスの速度分布vRvpvg φθ Gas: Ar, 400K

粒子速度依存性

散乱過程
実験室系

vp相対系

v p = v p (cos θ p cos ϕ p , cos θ p sin ϕ p , sin θ p )v p = v p (0,0,1)(vRを決定)v g = v g (cos θ g cos ϕ g ,cos θ g sin ϕ g ,sin θ g )

重心系

vR (0,0,1)→ vR (cos θ s cos ϕ s ,cos θ s sin ϕ s ,sin θ s )(ポテンシャル散乱)vR = v p − v g∞b / r2r01 − ( b / r ) 2 − U ( r ) / ( µv R 2 / w )θ s = π − 2∫dr高速粒子のガス中での拡散
Cu, B, La 原子 初期エネルギー: 5eV ガス:
Ar (5Pa, 400K)
 原子間ポテンシャル: Born-Mayer 型
U (r )=A exp(-br ) 追跡時間:
1~200 µs
 高速粒子:速度空間での時間発展
Cu/Ar 5Pa

T=2µs

T=10µs速度空間の時間発展

dN/dv [s/m]-2 Cu→Ar(400K)
10

10-3
10-4
10-50

T=50µs

5 µs
10 µs
20 µs
50 µs
100 µsMaxwellian

1000

Speed [m/s]2000速さ分布の時間発展（B, La）
10-2

dN/dv [s/m]B
10-35 µs
10 µs
20 µs
50 µs
Maxwellian

La

5 µs
10 µs
20 µs
50 µs
Maxwellian

10-4

10-502000

4000 0Speed [m/s]

10002000

拡散過程の取り扱い
 運動エネルギーが適当な閾値 (ガスの熱エネルギーなど)を下まわったスパッタ粒子は、「熱中性化(Thermalization)」したとみなす。
 熱中性化した粒子の拡散プロセスは Poisson 方程式 に支配される。
 粒子の熱中性化が起こった地点を Poisson 方程式での粒子の湧き出しとする。一つの粒子に対応する湧き出しは、３次元デルタ関数とする。
 境界での付着確率は 1 とする。

拡散過程(1) 拡散方程式(0 =)∂n p ( x)∂tρ (x) == D∆ n p ( x ) + ρ ( x )

1

NTH

TMCi3
δ
∑ (x − xi )境界条件: n p ( x B ) = 0

TMC: MC時間単位

単位面積への TMC あたりの到達数:∂ ∂
FTh = TMC D n p (x B )  ≡ n ⋅ grad ∂n ∂n拡散過程(2) MCとの整合η ≡ TMC Dn p を定義すれば、N TH

解くべき方程式:∆η ( x ) + ∑ δ ( x − x i ) = 0i

境界条件:

単位面積あたりの粒子到着数:

η (x b ) = 0∂ηdn Γ

両手法の計算結果の比較
Atomic Flux (a.u.)

5MC (10 )BEM

60402000

10203040

BE Node Numberまとめ: シミュレーションポリシー
容器、境界条件
 回転対称 ガスの温度、圧力は一様
 壁面での付着確率は 1
 粒子の飛び出し 位置
エロージョントラックの深さを反映
 Energy,角度Thompson の式、cosine 分布
 ガスとの散乱
 ガスの熱運動
Maxwell 分布
 散乱ポテンシャル Born-Mayer型 ( U(r)=Aexp(-br) ) 拡散過程 Poisson 方程式を境界要素法(一定要素モデル)で解く。
計算環境 IBM 互換 PC Debian GNU/Linux コンパイラ：GNU C & Fortran ライブラリ：LAPACK SLATEC Berkley Sparse Library ソースコード：

http://surf.ap.seikei.ac.jp/HISPUT/

シミュレーションの様子
Cu/Ar 5Pa

2µs

シミュレーションの様子
Cu/Ar 5Pa

20µs

シミュレーションの様子
Cu/Ar 5Pa200µs

粒子輸送過程の圧力依存性

膜厚分布の圧力変化
 膜厚分布→生産設備では重要
 銅のスパッタ製膜過程に対するシミュレーションの再現性を check する。
 Motohiro らの報告(1986)
圧力が上昇すると、基板ホルダ裏面への付着量が一旦増え、更に圧力を上げると減少する。
シミュレーションとの一致が良くなかったため、原因に関する考察も進んでいなかった。

元廣らの結果

T. Motohiro JVST A4 (1986) 189.

スパッタ製膜実験
DC マグネトロン
ターゲット： Cu
放電ガス : Ar
Erosion Track:x0.3 mm
0123

Radial Position (cm)

420 cmNormalized Flux (a.u.)

膜厚分布1.00.50-4Pressure: 2Pa
Calc (r=13cm)
Calc (r= 7cm)
Experiment-202

Pressure: 10Pa
Calc (r=13cm)
Calc (r= 7cm)
Experiment

4 -4-2

Substrate Position (cm)024

背面へのまわり込み
R/F Thickness Ratio

0.1040mm-35mm

0.05-30mm
01510Pressure (Pa)

熱中性化位置の圧力依存性

まとめ（1）: 粒子のまわり込み付着 低圧側 （熱中性化が起こらない）

圧力が上昇し、Ar ガスとの衝突が増えると、裏面へまわり込む銅原子が多くなる。
 高圧側 （熱中性化が dominant）

高圧になると熱中性化がよりターゲット近傍で起こるため、基板裏面までの距離が大きくなって到達率が下がる。高圧の極限では、裏面/表面の到達率の比が一定に近づく。

化合物スパッタ膜(LaB6)
 ターゲット組成と膜組成がずれる。ずれの大きさは圧力によって変化する。
 Rossnagel（1989）: 合金膜では、軽い元素が中間圧力で欠乏する傾向がある。
 LaB6
 構成元素がいずれも non-volatile ガスとの相対質量

M g M p /( M g + M p ) が Ar-B,Ar-La 間でほぼ等しい→エネルギー伝達が同程度

膜組成の圧力依存性
77B/La at.ratio(b) Simulation
665544Pos 5eV 10eV0- 510-1520-25(a) Experiment
30.5

151030.51Pressure (Pa)510

Scattering Angle (rad)

散乱角分布
33La

B
0.1eV

221eV

1 10eV

1000110eV 1eV 0.1eV2345 0Collision Parameter(Å)

12345Collision Parameter(Å)

B, La の拡散
Ar 2Pa

B

T=5µs

La

B, La の拡散
Ar 2PaB

T=50µsLa

B, La の拡散
Ar 2PaB

T=200µs

La

粒子のZ座標（2Pa）
Z position (cm)

86B

La

420

Ar 2 Pa-210-610-510-410-6Time (s)10-510-4

粒子のZ座標（5Pa）
Z position (cm)

86B

La

420Ar 5 Pa-210-610-510-410-6Time (s)10-510-4

まとめ(2): 組成変化
 低圧側 （熱中性化が起こらない）

ボロンの方がより後方に散乱されやすいため、圧力上昇とともにボロン組成が低下する。
 中間圧力（熱中性化が生じはじめる）

拡散過程には質量の差はあまり影響しないが、 La は前方に「打ち込まれて」から熱中性化して拡散を開始するため、Bよりも基板への到達確率が高い。
 高圧側 （熱中性化が dominant）

高圧になるほど熱中性化がターゲット近傍で起こり、拡散開始地点の差が小さくなって、組成は回復する。

スパッタ原子の空間密度
ターゲット
スパッタリング

プラズマ
粒子輸送過程励起・イオン化過程
デポジション(フラックス・エネルギー・角度）

膜物性

 膜物性の改善のために、付着粒子のエネルギー・角度を制御
 “ionized sputtering”
プラズマ中での励起・イオン化
 粒子の空間密度・滞留時間などの評価が必要

MC 法と粒子密度
 空間にセルを配置、それぞれを粒子が通過する時間を求める。
 密度はセルを出入りする粒子のフラックスと滞留時間の積に比例。

通過時間Tp密度の導出Simulation

Experimenty 壁面付着数

y 成膜速度y 試行粒子数y 実時間 Tsy 領域通過時間 Tp

y 粒子数 N→密度 n

Substrate

4cmTarget

DC マグネトロンスパッタ
Target Cu / Gas Ar
Power 100W

Deposition Rate (nm/s)

成膜実験
0.2

0.1

01

510

Pressure (Pa)

計算結果

最近の実験報告

M. F. Dony, et al., JVST A18(3) pp.809 (2000)
Al 金属やガラスを Ar でスパッタ、吸収分光で密度を評価。
3.5cmφの Al ターゲットに対し、RF 30〜40W、Ar 圧力
40〜80 mTorr で 3-4×1011cm-3 程度。

 Y. Andrew, et al., JAP 88(6) pp.3208 (2000)

Cu/Ar イオン化スパッタの吸収分光。 15cmφのターゲッ
トに DC 0~200W で、10〜30mTorr で 1012~1013cm-3。
イオン化用の RF コイル電力を上げると、Cu 密度も上昇。

基底原子による発光の自己吸収
Cu I

3.817eV

510.6nm
324.8nm1.389eV
0eV

n (m-3)λ 324(m)

1016
1017
1018

1.4×10-1
1.4×10-2
1.4×10-3

 Cu I 324.8 nm の発光は基底状態への遷移なので、周囲に Cu原子があると吸収される。
 Cu 原子密度が大きくなると、同じ始状態の510.6 nm の発光が相対的に大きくなる

Cu/ArプラズマでのOES2 Pa

Intensity (counts)

Intensity (counts)

2000
1000
02000

20 Pa

1000

10 3

10 2324.8nm
510.6nm

10100 W
0200300400

Wavelength (nm)500111020 Pa
1050 100

Pressure (Pa) Power (W)まとめ(3) 空間密度分布
 圧力が高くなるにつれて密度が上昇した。

平均自由行程の減少と、速やかな減速とによって拡散係数が小さくなるため。
 Cu 原子からの発光の圧力・電力依存性計算によって求めた粒子密度と、発光線の自己吸収効果によってうまく説明できた。

まとめ スパッタ製膜プロセスのモンテカルロシミュレー

ションを拡張し、数Pa以上の高い圧力でもシミュレートが可能になった。
 ガス熱運動のMCシミュレーションへの取込み
 拡散輸送のPoisson方程式による取扱い

 膜厚分布や化合物膜の組成変化の圧力依存性

は、熱中性化前後のballisticな輸送とdiffusiveな輸送の競合で説明できる。
 空間密度の導出など、スパッタ粒子の空間挙動に対してもMC法が利用できそう。