分子動力学法による格子間シリコン原子の拡散パス

90年代半ばまでは、我々にできることは、測定可能な“デバイスの電気特性”と“一次元不純物分布”だけ から、裏に隠れた原子レベルの振る舞いを「想像して」モデル化することであって、その想像の真偽・精度 を直接実験により確認することは困難であった。ところが、計算機の進歩と物理化学的計算手法の進歩に伴 い、密度汎関数理論に基づく第一原理分子動力学計算手法が登場し、この「想像」を計算機内であたかも起 こっているように調べることが可能になってきている。

一方、「設計の効率化」を目的としたモデリングは、なるべく計算効率の良いシミュレーションツールに 帰着すべきであり、その意味では分子動力学法は直接の設計ツールの基礎を得る手法と位置づけられるべき ものである。本節では、シリコン中のイオン注入と不純物拡散についての原子レベルのプロセスシミュレー ションの基礎的データを提供する分子動力学計算についてその有効性を議論する。

3.6.1 古典的分子動力学 (MD) 計算

格子間シリコン原子拡散

1985年StillingerとWeberによりシリコン結晶の古典的な原子間ポテンシャルが発表され［110］、シリコ ン結晶の熱・機械的性質が計算機上で再現されるようになった。この古典的ポテンシャルを用いてシリコン 中の点欠陥の拡散の解析も行われるようになっているが［111］［112］［113］、本研究では、より詳細な、イオ ン注入工程やシリコン中の点欠陥の拡散の計算のため、この分子動力学計算の方法を以下に述べる方法で改 良し、原子数にほぼ比例する程度まで計算時間を短縮させ、より実用的な計算を可能にした［113］。

Stillinger-Weberポテンシャル（以下SWポテンシャル）に基づく分子動力学計算とは、解析的に表わされ

た２体・３体原子間ポテンシャルに基づくニュートンの運動方程式を数値積分する計算である。

mid²ri

dt² =Fi=−∇i

XN

j6=i

Φij (3.58)

ここで,m、rはそれぞれ各原子の質量と位置座標、Nは原子数を表わす。SWポテンシャルは次のような関 数形で表わされる。

Φij =X

i<j

φ2(rij) + X

i<j<k

φ3(ri,rj,rk) (3.59)

φ2(rij) =εf(rij/σ) (3.60)

φ3(ri,rj,rk) =εg(ri/σ,rj/σ,rk/σ) (3.61) ここで,

f(x) =

½A(B/x^p−1/x^q) exp[1/(x−a)] x < a

0 x > a (3.62)

g(xi,xj,xk) =h(xij, xik, θjik) +h(xji, xjk, θijk) +h(xki, xkj, θikj)

h(xij, xik, θjik) =









λexp(γ/(xij−a)+

γ/(xik−a)(cosθjik+ 1/3)²

· · ·xij < aandxjk< a 0 · · ·xij > aorxjk> a

(3.63)

上記のパラメータとしてA= 7049556277,B = 0.6022245584,p= 4,q= 0,a= 1.8,λ= 21.0,γ= 1.2, ε= 3.4723×10⁻¹⁹J,およびσ= 0.20951nmがシリコンの格子定数、融点、体積弾性率等を良く再現する セットとして文献［110］に挙げられている。

SWポテンシャルによる分子動力学計算プログラムを組む際、原子３個の３体の組み合わせを探索する計 算量は、全原子数をN とすると、総当り探索では_NC₃と書ける。この場合、計算の規模を大きくする、すな わちNを増加させると計算時間はN³(∼ O(N³))オーダーで増加することになる。しかしながら、式3.63 中のh(xij, xik, θjik)の項はxij,xik> a, (a:カットオフ長)の場合のみゼロをとらない。よって、カットオ フ長内の隣接原子についてhの構成要素として各原子に関しseparable form［114］を用いたデータテーブル 化が全体の計算を効率化する。これに必要な計算量は全体から２体の組み合わせを探す分NC2に減る。３ 体のポテンシャルφ3は、カットオフ長a内の原子を探すO(N²)のオーダーの手間で、各原子のseparable formから合成することができる。

更に隣接原子の探索を、着目原子の隣接ブロック内でのみ探索するブロックサーチのアルゴリズムにより 行った。最終的にSWポテンシャルを組み上げる計算量はO(N)のオーダーまで減らした。すなわち、あら かじめ領域全体をカットオフ長よりもやや大きい直方体でブロック分けを行ってからブロック内の原子数Nb

のみで探索を行う計算量は、N/Nbに比例する。ブロックサーチは計算したい系が大きくなるほど計算を効 率化できる。Separable-formによる計算順序変更による効率化とブロックサーチによれば、最終的にSWポ テンシャルを組み上げる計算量はO(N)のオーダーまで減らせる。典型的な計算時間の例としては、NEC

製EWS4800(R10000)で、500原子ユニットセルにおける1タイムステップの計算時間は、総当りサーチで

は12.0秒を要するが、本研究の場合、1.0秒以下であった。

図3.44にシリコン中の格子間シリコンの拡散をSWポテンシャルMDで計算した例を示す［113］。

!"#

#

$! %&

'

()+*-,+,-.0/213+46587 9:)+*-,+,-.8/213+46587

;+)-*+,-,+.</213=46567

図3.44:シリコン結晶中の格子間シリコン原子(self-interstitial)の1073Kにおけるマイグレーションの 2nsまでの計算結果

分子動力学計算のアンサンブルとしては、NVE (粒子数N,体積V,エネルギーE,がそれぞれ一定）とNVT

（温度T一定）を用意し、時間ステップの更新には、４次のギア法を用いた。 格子間シリコン原子の拡散の 分子動力学計算には、シリコン結晶格子5×5×5のシリコン原子からなるユニットセル（原子数1001）を

用い、イオン注入の分子動力学計算は14×14×14格子分のユニットセルを用いた。例えば、シリコン結晶 について、1000原子程度の単位セルで数ナノ秒の熱拡散時間の計算が、単体ワークステーション上で数百時 間程度の計算時間で可能である［113］。

熱拡散の分子動力学計算結果から次式のアインシュタインの関係式で拡散係数を抽出することができる。

D_I =X

(x−x₀)²/6t (3.64)

ここで、x0は原子の初期位置、xは時間t後の全変位量である。異なる一定温度(NVTアンサンブル)の計 算を少なくとも３条件行ってアレニウス型の表式も得ることができる。例えば図3.44の例では、式(3.65)の ような結果となる。

DI = 0.039 exp(−0.89/kT)[cm²/sec] (3.65) ここで得られた活性化エネルギー(0.89eV)は、異なる計算条件の他の報告値(0.94 eV［112］または0.9 eV

［111］)ともほぼ一致している。ただし、有限な温度下での大規模なMD計算は起りうる全ての物理的過程 を含んだ模擬的な結果であって、例えばこの結果を基に詳細な格子間シリコン原子の拡散過程を理解したと は言い難い。しかしながら、このMD計算から得られた活性化エネルギー0.89eVが、ミクロ的にどのよう な過程から生じているのかは、以下に述べる反応経路の同定という手法を使って調べることが可能である。

次節でそれについて述べる。

格子間シリコン原子拡散パスの同定

図3.44からは、格子間シリコン原子の拡散がある安定サイトを渡り歩くホッピングのような振る舞いであ ることが判る。格子間シリコン原子１個を含む有限温度下の熱拡散の分子動力学計算結果から、各時間刻み 毎の格子間シリコン原子周りの構造から出発して絶対温度零度まで構造緩和させていくと、ある典型的な格 子間シリコン原子の安定原子配置が得られる。拡散は異なる安定構造サイト間をホッピングしていく形で起 り、ホッピング間にいわばエネルギーの鞍点(saddle point)の構造がある。

Nastarらは、SWポテンシャルを用いて格子間シリコン原子の安定サイトとしてh110iダンベル構造が

h110i方向に遷移（移動）するものと仮定して、独自の鞍点探索アルゴリズムを用いて遷移エネルギー1.62eV を得たと報告している［115］が、明らかに前述のMD計算から得られる活性化エネルギー0.89〜0.94eVより 高すぎる。つまり、Nastarらの「仮定」した反応パスは、(SWポテンシャルが語る)格子間シリコン原子の 拡散の素過程では無い。

鞍点構造を見つけるアルゴリズムとして制限付き構造緩和手法が試された［113］。参考文献［116］には着目 したある一つの原子に働く力を式(3.66)のように制限することで、鞍点のエネルギーを得る方法が提案され ているが、格子間シリコン原子の遷移には周りの複数原子の関与が大きいので、この方法を全原子に拡張し て用いた。

F∼=F−(F,4R)4R

|4R|² (3.66)

ここで、4Rは終状態構造への微小強制変位ベクトル、Fは各原子に働く力で、F^∼は制限された力である。

あらかじめ行った分子動力学計算から、拡散パスは大まかには、h110iダンベル構造がh110i方向よりも ねじれた方向のサイトへ動くことが見て取れるので、この始状態と終状態の構造について、制限付き構造緩 和法により求めた反応経路に沿った全エネルギー変化を図3.45に示す［113］。

0.0 0.2 0.4 0.6 0.8 1.0 1.2 Pathway length (nm)

Energy (eV)

3.0 3.5 4.0 4.5

5.0

Extended <110>

dumbbell Extended <110>

dumbbell

Em (a)

(b) (c) (d)

(e)

(f)

図3.45: SWポテンシャルによる格子間シリコン原子拡散パス上の全エネルギー変化

図3.45から読み取れる遷移エネルギーバリアEmの値は、0.9eVであって、分子動力学計算から抽出され る拡散係数の活性化エネルギーに相当するものとなっている。このように、有限温度の分子動力学計算と安 定構造と鞍点構造計算を組み合わせて、諸現象の「モデル」を素過程として組み立てることができる。

3.6.2 半経験的量子力学的ＭＤ法

Stillinger-Weberポテンシャル以外にも、古典的な分子動力学用の原子間ポテンシャルにはTersoffポテン

シャルなど様々な提案がありシリコンの非晶質化等が調べられている［117］。しかし、SWポテンシャルが 良く再現する融点の実験結果をTersoffポテンシャルは再現しない。一方、Tersoffポテンシャルは、点欠陥 の局所構造は、SWポテンシャルよりも第一原理計算の結果を良く再現する。このように機械的な構造のみ をフィッティングしたポテンシャルによる古典的なMD法には、その構造の電子状態はどうなっているのか という本質的な疑問がどうしても付きまとう。原子間力を量子力学に基づく電子状態を解いて求める必要が あるが、かといって第一原理計算では現実的な計算規模には到底収まらない。そのような場合、ある限られ た特徴的な第一原理計算結果とうまくフィッティングされ、かつ比較的簡便に行える電子状態の計算モデル

としてTight-bindingモデルがある。これを用いれば、実用的な問題に対しても半経験的な量子力学的MD

解析が可能になる。

Di-interstitialの拡散

イオン注入後の拡散で問題になる過渡的な増速拡散を引き起こす過剰格子間シリコン原子の寿命を律速し ているのは、格子間シリコン原子の凝集形態である{311}欠陥であると言われている［76］。{311}欠陥の成 長機構も分子動力学的手法で議論されている［118］。また{311}欠陥の前駆体として、格子間シリコン原子 ２個のペアであるdi-interstitialが着目されている［119］。Stillinger-Weberポテンシャルによる分子動力学計