修士論文

半整合界面の構造安定性ならびにミスフィット転位に関する分子動力学的研究

指導教員：屋代 如月

鈴木 雄風2009 年 2 月
神戸大学大学院 工学研究科 博士課程前期課程 機械工学専攻

Molecular Dynamics Study on Stabilityof Semi-coherent Interfaces and Misﬁt Dislocations

February 2009
Department of Mechanical Engineering,
Graduate School of Engineering,
Kobe University, Kobe, Japan
Yuuhu SUZUKI

要約
本研究では，異種材料界面に形成される転位をテーマに，Ni 基単結晶超合金 γ/γ ′ 界面
ならびに半導体デバイスの Si/GaAs 界面について，種々の分子動力学シミュレーショ
ンを行った．まず，Ni 基超合金の γ/γ ′ 界面上のミスフィット転位と，γ 相を進行してきた刃状およびらせん転位の切り合い，ならびに，外力下でのミスフィット転位の運動について検討した．ミスフィット転位は，異符号側に余剰原子面を持つ刃状転位や，同一のバーガースベクトルを持つらせん転位とは転位論により説明可能な挙動 (バーガースベクトルの可算則) を示した．しかしながら，ミスフィット転位と同符号側の刃状転位が十字に切り合うとき，ミスフィット転位を陵線とする三角屋根状の 2 つのすべり面上に新たな転位を生じ，衝突した刃状転位と三角屋根の片方に発生した転位の交点が，それぞれのすべり面が交差する線 (V 字の溝) 上をスライドするメカニズムを見出した．また界面に平行にせん断を与えるシミュレーションを行い，[100] 方向へのせん断では，界面転位網が界面上に多数の点欠陥を残しながら全体的にすべること，[110] 方向へのせん断では，転位網の横糸部分のみ界面をすべること，などを示した．
Si/GaAs 界面の検討では，Si と GaAs それぞれに Tersoﬀ 型ポテンシャルが提案されているものの，界面での取り扱いが不明であることから，まず単位格子を接合したスーパーセルによる解析を行い，界面エネルギーを第一原理計算により検証した．その後，

Si 基板上に GaAs ハットクラスターを設置し，構造緩和するシミュレーションを行い，ポテンシャルパラメータの取り扱いによる界面欠陥構造の違いについて議論した．Siの表面構造を再現する Tersoﬀ 型ポテンシャル (T2) を用いると，クラスターの縁または下の界面に Si の再構成による欠陥を生じるが，弾性係数を良好に再現する T3 を用いると，Si，GaAs それぞれは単位結晶でのダイヤモンド構造から大きく変化せず，界面上下の格子を比較したときにのみズレとしての転位が存在する．最後に，無限平面とした Si/GaAs 界面に，界面法線方向へ引張を与えるシミュレーションを行った．Si にT2 を用いた場合は，Si 基板上の初期欠陥を起点として，界面から {111} 面に沿ったき裂が生じて破断した．一方，T3 を用いた場合は，界面から離れた GaAs 層内に {111}
面のへき開を生じた．

Summary
The main ofjective of this study is the atomic-scale detail of dislocations at theinterface between diﬀerent materials. From such point of view, the γ/γ ′ interface in
Ni-based superalloys, the Si/GaAs interface of semiconductor device are investigatedby molecular dynamics simulations. First we have observed the interaction between
misﬁt dislocations on the γ/γ ′ interface and edge/screw dislocations gliding in a monolattice, and motion of misﬁt dislocation network under the shear force. It is revealed
that the reaction can be basically explained by the dislocation theory, i.e. Burgers
vector problem. However, we have also found new mechanisms in the cross cuttig by
an edge dislocation which has extra plane on the same side as the misﬁt dislocation,i.e. same sign. New dislocations are generated on the delta-roof-like two slip planes
from the peak line of [110] misﬁt dislocation. Furthermore, the edge dislocation and
one dislocation on the delta-roof have nodal point gliding on the V-gutter between two
slip planes. We have also applied the [100] and [110] shear on the interface. Under
the [100] shear, the whole network glides dragging many defects. On the other hand,
only the weft of the network glides without leaving defects, under the [110] shear.
The critical stress for the constant slip of the interface is slightly diﬀerent for the[100] and [110] shear, due to the anisotropy of the lattices. For the atomic simulationof Si/GaAs interface, we have an ambiguity for treatment of potential parameters
between Si/GaAs, although there is a famous potential function for Si and for GaAs
by Jeﬀerson. Thus we have ﬁrst checked the Si/GaAs interfacial energy by DFT
ab–initio calculation (VASP) with a super cell stacking the Si and GaAs unit cells.
Then, we have performed molecular dynamics simulation on GaAs hut cluster on Si
substrate and discussed the defect structure at the interface. If we use the Tersoﬀ
potential for surface reconstruction of Si (T2), we can ﬁnd dissorders of Si substrateat the edge and interface under GaAs cluster. The Tersoﬀ parameter for the proper
elastic coeﬃcient (T3) leads no disorder on the interface and the crystal lattices ofSi and GaAs independently remain at the initial diamond structure. We can ﬁndlattice mismatch only if we compare the upper and lower crystal lattices. Finally wehave applied tensile loading on an inﬁnite Si/GaAs interface. We have observed crack
initiation and propotion in GaAs phase from the initial defect on Si substrate with theT2 parameter. On the other hand, we have observed cleavage cracking in the middleof GaAs layer, with T3 parameter.

目次
第1章 緒 論1

Ni 基超合金中の γ/γ ′ 界面 . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2 半導体デバイス中の Si/GaAs 界面 . . . . . . . . . . . . . . . . . . . . .

2

1.3 本論文の構成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.1

第 2 章 分子動力学法

52.1 分子動力学法の概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

2.2 原子間ポテンシャル . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2.3 原子埋め込み法ポテンシャル. . . . . . . . . . . . . . . . . . . . . . .

7

Tersoﬀ 型ポテンシャル . . . . . . . . . . . . . . . . . . . . . . . . . . .

92.5 速度スケーリング法 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

2.6 高速化手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

2.4

第 3 章 Ni 基超合金 γ/γ ′ 界面における転位挙動

17

3.1 シミュレーション条件 . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

3.1.1

解析モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

3.1.2

刃状転位シミュレーション . . . . . . . . . . . . . . . . . . . . .

18

3.1.3

らせん転位シミュレーション . . . . . . . . . . . . . . . . . . .

20

3.1.4

せん断シミュレーション . . . . . . . . . . . . . . . . . . . . . .

21

3.2 シミュレーション結果および考察 . . . . . . . . . . . . . . . . . . . . .

23

3.2.1

刃状転位シミュレーション 1 . . . . . . . . . . . . . . . . . . . .

23

3.2.2

刃状転位シミュレーション 2 . . . . . . . . . . . . . . . . . . . .

28

3.2.3

らせん転位シミュレーション . . . . . . . . . . . . . . . . . . .

33

3.2.4

せん断シミュレーション . . . . . . . . . . . . . . . . . . . . . .

36

3.3 結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

i

目次
第 4 章 Si/GaAs 界面における転位挙動

ii
43

4.1 界面エネルギーの評価 . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

4.1.1

三軸等方膨張・収縮を受ける Si・GaAs 単相のエネルギー評価 .

43

4.1.2

界面におけるポテンシャルパラメータの扱い . . . . . . . . . . .

47

4.1.3

Tersoﬀ 型ポテンシャルによる界面エネルギーの評価 . . . . . . .

50

4.2 界面原子構造の評価 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

4.2.1

シミュレーション条件 . . . . . . . . . . . . . . . . . . . . . . .

54

4.2.2

シミュレーション結果および考察 . . . . . . . . . . . . . . . . .

56

4.3 引張シミュレーション . . . . . . . . . . . . . . . . . . . . . . . . . . .

62

4.3.1

シミュレーション条件 . . . . . . . . . . . . . . . . . . . . . . .

62

4.3.2

シミュレーション結果および考察 . . . . . . . . . . . . . . . . .

64

4.4 結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67

第 5 章 結論

69

参考文献

71

関連発表論文・講演論文

76

謝辞

98

第1章
緒 論

1.1

Ni 基超合金中の γ/γ ′ 界面

Ni 基単結晶超合金は大型ガスタービンやジェットエンジン等に広範に用いられる主要な耐熱材料である．その内部構造は，Ni3 Al を主成分とする γ ′ 相を，Ni 固溶体である γ 母相中に立方体形状で格子状に微細析出させた特徴的な形態を有する．γ および

γ ′ 相の結晶格子は非常に近い関係にあり，格子ミスフィットが非常に小さいため，γ とγ ′ 相は整合界面を形成する．このため，このような内部構造を有しながらも全体としては単結晶である．一方，クリープ条件下では，γ ′ 相が引張軸に対し平行または垂直方向に粗大化する (ラフト化)(1) ．ラフト化すると，γ/γ ′ 界面の面積が初期構造に比べて極めて大きくなるため，弾性ひずみ場を緩和するように規則正しいミスフィット転Fig.1.1 Typical morphology of the microstructure of Ni–based superalloy11.2 半導体デバイス中の Si/GaAs 界面

2位が界面に生じる．このミスフィット転位の間隔がクリープ寿命と密接に関連しているとの報告もあり(2)–(4) ，第 5 世代の超合金では転位網間隔を積極的に制御することで耐熱性能を向上させるという試みもある(5) ．
顕微鏡技術の進歩により，転位の挙動の動的な観察が可能となっているが，薄膜状試料に限られるなど限界があり，数値シミュレーションによる解明が期待されている．
多数の転位の集団的挙動を扱う手法として離散転位動力学(6), (7) があるが，転位論に基づいて定式化されているため，結晶の格子構造の差を考慮できるものではない．半整合界面上のミスフィット転位は，界面上下の原子構造を比較してはじめて定義される．
すなわち，界面上下の相を見た場合，原子面が不足または余っているわけではない．ミスフィット転位の応力場等は，転位論により評価できるが，転位芯同士の相互作用 (切り合いなど) については，原子レベルからの検討が必須である．そこで，分子動力学による検討が挙げられる．Ni 基超合金に関する研究には Ni/Ni3 Al の整合界面および Al
濃度を傾斜させた界面の構造安定性および弾性特性に関する検討(8) ，整合界面近傍での転位挙動(9), (10) などの他，半整合界面上のミスフィット転位については，様々な結晶
方位で接合した Ni/Ni3 Al 界面に形成される転位網の構造およびエネルギーに関する検討(11) ，引張または圧縮による構造変化(12) ，さらに γ 相を進行するプリズマティック転位ループとミスフィット転位の相互作用(13) 等の研究がある．しかしながら，γ 相を進行する 1 本の転位とミスフィット転位の相互作用を詳細に議論したものはない．

1.2

半導体デバイス中の Si/GaAs 界面

直接遷移型である III–V 族化合物半導体は，半導体材料中を移動する電子移動度がシ
リコンの数倍〜数十倍もあるため非常に注目されている．しかし，化合物半導体は脆く大型化が困難なため，良質な Si 基板上にエピタキシャル成長させて用いられている．

Si 基板上に GaAs を成長させると，GaAs ハットクラスターが形成される(14) ．はじめ
は僅かなひずみエネルギーがあるが，転位が入るほどエネルギーは高くならず，Si とGaAs が整合する(14)–(17) ．しかし，GaAs をさらに成長させるとひずみエネルギーが大1.3 本論文の構成

3

きくなり，整合するよりも転位を導入した方がエネルギーが低くなるため，はじめは長方形界面の片側の面に転位が入り，その後，界面上に周期性を持ったミスフィット転位が導入される(18) ．ミスフィット転位は界面に形成するが，材料中には多数の欠陥が存在する．例えば，GaAs 層へ伝播する欠陥として双晶境界や積層欠陥がある(15) ．この積層欠陥は結晶成長時から存在する貫通らせん転位またはミスフィット転位からの拡張転位と，熱処理時の外部応力によって引き起こされた滑り転位の 2 種類が存在す
る(19) ．長寿命な優れたデバイスを作成するためには，GaAs 層をエピタキシャル成長
させる際に，転位の起点となる界面の欠陥を制御することが重要である．

Si と GaAs の格子定数の差は 4%もあるため，ミスフィット転位の導入を回避することは困難である．しかし，ミスフィット転位から伝播する転位等の欠陥を限りなく減らすことにより，長寿命の優れたデバイスの創成が期待できる．その基礎として，ミスフィット転位を制御すること，およびミスフィット転位の生成メカニズムを解明することが重要となる．Si に関する数値シミュレーションは Ni 基超合金のそれに比べて非常に多くなされている(20)–(23) ．エピタキシャル成長過程に関する研究には，基板温度や原子の初期入力位置の違いによる応力評価及びシリコン薄膜の真性応力に関する検討(20) ，InAs/GaAs(110) ヘテロ界面の芯構造と電子状態及びミスフィット転位の構造
が自己組織化量子ドット形成に与える影響に関する検討(21) ，GaN(0001) 表面で吸着窒
素原子の過多条件による結晶成長の促進に関する検討(22) ，GaAs 表面での吸着とエネ
ルギー的に不利な Ga あるいは Si ダングリングボンド中の電子数の関係(23) などの研究
がある．

1.3

本論文の構成本研究では，Ni 基超合金 γ/γ ′ 界面上のミスフィット転位および Si/GaAs 界面上のミ
スフィット転位について，原子レベルから知見を得ることを目的として種々の分子動力学シミュレーションを行った．以下，各章の概略を示す．
第 2 章では解析手法について述べる．はじめに，分子動力学法の基礎方程式を示し，

1.3 本論文の構成

4

本研究で原子間相互作用の評価に用いた原子埋め込み法ポテンシャル及び Tersoﬀ 型ポ
テンシャルの概要について説明する．また，大規模分子動力学計算の手法である領域
分割による高速化について説明する．第 3 章では，γ/γ ′ 界面上のミスフィット転位と，刃状転位またはらせん転位との転位芯同士の相互作用を分子動力学シミュレーション
により詳細に検討する．また，同モデルを用いてせん断シミュレーションを行い，その時のミスフィット転位の挙動についても検討する．第 4 章では，はじめに Si/GaAs
界面におけるエネルギーの評価を行う．次に Si 基板上に GaAs を乗せたモデルを用いて界面上に生成する欠陥原子について検討する．最後に，第 5 章で本研究の総括を述
べる．

第2章
分子動力学法2.1

分子動力学法の概要

分子動力学法 (Molecular Dynamics Method; MD) は，系を構成する個々の原子についてニュートンの運動方程式

mid2 r i
= Fi
dt2

(2.1)

を作成し，これを数値積分することによって全原子の運動を追跡する手法である．ここで t は時間，r i ，mi はそれぞれ原子 i の位置ベクトル及び質量である．原子 i に作用する力 F i は系全体のポテンシャルエネルギー Etot の空間座標についての勾配ベクトルから次式のように求められる．

Fi = −

∂Etot
∂r i

(2.2)

式 (2.1) の数値積分には，Verlet の方法が簡便で高精度が得られるため MD 法ではよく用いられる．時刻 t±∆t での原子 i の座標 r i (t ± ∆t) を Taylor 展開すると

r i (t ± ∆t) = r i (t) ± ∆t

dr i (∆t)2 d2 r i (∆t)3 d3 r i+±+ O((∆t)4 )dt
2 dt2
3! dt3(2.3)となる．両式の和をとり式 (2.1) を代入すると

r i (t + ∆t) = 2r i (t) − r i (t − ∆t) + (∆t)2

F i (t)+ O((∆t)4 )
m

(2.4)

を得る．これより，時刻 t − ∆t と t における全原子の位置を既知として，時刻 t + ∆t
における任意の原子 i の位置を求めることができる．

5

2.2 原子間ポテンシャル

2.26

原子間ポテンシャル

系のポテンシャルエネルギー Etot の評価は，以下の 3 つに大別される．

(1) 経験的ポテンシャル(2) 半経験的ポテンシャル(3) 非経験的手法 (第一原理計算)
経験的ポテンシャルは，量子力学の厳密な理論に基づいて決定されるのではなく，ポテンシャルを微分可能な未定係数を含む簡単な関数形で仮定し，従来の実験的事実に合致するようにその未定係数が決められる．半経験的ポテンシャルは，密度汎関数論
より導出される形で定義されるが，そのポテンシャルパラメータは平衡状態でのマク
ロな物性値や，あるいは ab-initio な計算により求められた値に対してフィッティング
される．非経験的手法とは，従来の特性値などを一切用いず，原子核の位置ならびに種類のみを必要情報とし，各時刻における電子状態を量子力学に基づいて解くことで，遂次原子に働く力を精密に評価する手法である．

MD 法において，原子 i に作用する力 F i は系のエネルギー Etot の空間微分によって求めるため (式 (2.2))，系のポテンシャルエネルギー Etot をいかに精度よく評価するかが重要となる．(3) の第一原理分子動力学法は，計算量が極めて膨大になるため，変形・破壊のような多数の原子の動的挙動への直接的な適用は困難である．そこで，原子間相互作用を簡略評価する (1)，(2) の原子間ポテンシャルが通常用いられる．本解
析の原子間ポテンシャルでは，Ni 基超合金には Daw，Baskes らによって提案された原子埋め込み法 (Embedded Atom Method; EAM)(24), (25) ，Si/GaAs には Tersoﬀ 型ポテ
ンシャル(26)–(28) を用いている．

2.3 原子埋め込み法ポテンシャル

2.37

原子埋め込み法ポテンシャル

EAM は金属中の多体効果を良好に再現することから広く用いられている．EAM では密度汎関数理論に基づき，まず金属材料における系のポテンシャルエネルギー Etot
は原子を価電子雲中に埋め込むエネルギーと原子間の 2 体間相互作用の和で与えられるとする．さらに，埋め込みエネルギーは埋め込む位置の電子密度にのみ依存する (局所密度近似) と仮定することによって，系全体のエネルギーは次式のように表わされる．

Etot =N∑

F (ρ̄i ) +

i

N
N ∑
1∑
ϕ(rij )
2 i j(̸=i)(2.5)

ここで，ρ̄i は原子 i の位置における電子密度，F (ρ̄i ) は電子密度 ρ̄i の位置に原子を埋め込むエネルギー，ϕ(rij ) は距離 rij 離れた原子 i と j のクローン相互作用である．また，ρ̄i は周囲の原子 j からの寄与 ρ(rij ) の重ね合わせで与えられると仮定しρ̄i =neighbor∑ρ(rij )

(2.6)

j(̸=i)

で評価する．

2 種類以上の原子を含む系では，系のエネルギーは次式で表される．
N
∑

N ∑
N
1∑
Etot =
Fti (ρ̄i ) +ϕt t (rij )2 i j(̸=i) i j
i

ρ̄i =

neighbor∑ρtj (rij )

(2.7)

(2.8)

j(̸=i)

ここで Fti (ρ̄) は原子種 ti の埋め込みエネルギー関数，ϕtitj (r) は原子種 ti と tj 間の 2 体間相互作用，ρtj (r) は原子種 tj の電子密度関数である．

Ni, Al の 2 原子系への適用を考える場合，FNi (ρ̄)，FAl (ρ̄)，ρNi (r)，ρAl (r)，ϕNiNi (r)，ϕAlAl (r)，ϕNiAl (r)，の関数形が必要となるが，′
ϕNiNi(r) = ϕNiNi (r) − 2 gNi ρNi (r)(2.9)

FNi′ (ρ̄) = FNi (ρ̄) + gNi ρ̄Ni

(2.10)

2.3 原子埋め込み法ポテンシャル

8なるパラメーター gNi (または gAl ) を導入して関数形の変換を行っても，式 (2.7) で Ni,

Al それぞれの純金属系におけるエネルギーは変化しない．さらに′(r) = sAl ρAl (r)ρAl(2.11)

FAl′ (ρ̄) = FAl (ρ̄/sAl )(2.12)のように Al の価電子密度をスケーリングするパラメーター sAl を導入して Ni 単原子系の電子密度の範囲に合わせる変換を施しても不変である．したがって，Ni, Al それぞれ単原子系で FNi (ρ̄)，FAl (ρ̄)，ρNi (r)，ρAl (r)，ϕNiNi (r)，ϕAlAl (r) を決定した後に，Ni–Al
合金系のエネルギーに対し ϕNiAl (r)，gNi ，gAl ，sAl を最適化することにより，Ni–Al 合
金系のエネルギーを正確に表すポテンシャル関数が決定できる．
本解析で対象とする Ni ならびに Ni3 Al については，Voter らが Ni，Al，Ni3 Al それ
ぞれ単結晶に対して昇華エネルギー，空孔形成エネルギー，弾性定数，格子定数等へのフィッティングを行い，ρ(r)，ϕ(r) の関数形として，

ρ(r) = s r6 (e−β r + 29 e−2 β r )(2.13)

ϕ(r) = D {1 − exp[−α (r − R)]}2 − D − 2 g ρ(r)(2.14)

を提案している(29)–(31) ．式中のパラメーターの値は表 2.1，2.2 に示した．

Table 2.1 Potential parameters for ρ(r)NiAl

β (Å−1 )3.64083.3232s
1.00000.6172

Table 2.2 Potential parameters for ϕ(r)
Ni-Ni
Al-Al
Ni-Al

D (eV)
1.55353.77603.0322α (Å−1 )1.77281.48591.6277R (Å)
2.2053
2.11762.0896g (eVÅ3 )6.5145-0.220502.4 Tersoﬀ 型ポテンシャル

9

埋め込みエネルギー関数 FNi (ρ̄)，FAl (ρ̄) については，原論文中でその具体的な関数形は示されていないが，Foiles(32) が提案している Rose らの純金属系の凝集エネルギー関数(33) を用いる方法で数値的に求めることができる．本研究では，この方法により埋め込みエネルギー関数を 3 次のスプライン関数でフィッティングした．フィッティン
グ範囲は 0.0 ≤ ρ̄ ≤ 1.0 で，スプラインノードの間隔は ∆ρ̄ = 0.01 とした．

2.4

Tersoﬀ 型ポテンシャル

Si や GaAs 等の共有結合材料は sp3 の強い方向性を持った結合を有するため，単純な 2 体間ポテンシャルでは表現できない．そこで結合角の効果を取り入れた多対ポテ
ンシャルを用いる必要がある．共有結合材料の代表的なポテンシャルである Tersoﬀ 型
ポテンシャル(26)–(28) は，経験的ボンドオーダーポテンシャルと呼ばれ，配位数や結合角などの周りの環境に合わせて結合状態を変化させている．また，Si のポテンシャル
パラメータには，多形態のシリコンの凝集エネルギー等がフィッティングされており，様々な条件下での原子構造を柔軟に表現できるようにしている．
本解析の MD シミュレーションで用いる Si および GaAs は，共に共有結合材料であるため，Tersoﬀ 型ポテンシャルを用いた．Tersoﬀ 型ポテンシャルでは，系全体の結合エネルギーを各原子ごとのエネルギーの和で表わされる．

Etot =

N
∑

Ei =

i

N ∑
N
1∑
Eij
2 i j(̸=i)(2.15)

各原子ごとのエネルギー Eij は次のように表わされる．
[

Eij = fc (rij ) fR (rij ) + b∗ij · fA (rij )](2.16)

fR = A exp(−λ1 rij )(2.17)

fA = −B exp(−λ2 rij )(2.18)

1b∗ij = (bij + bji )2(2.19)

2.4 Tersoﬀ 型ポテンシャル


1,

 [


fc (rij ) =






1 − sin(

π (rij − R)2D

10

: rij < R − D)]: R − D < rij < R + D/20,(2.20): rij > R + D

ここで fR (rij ) は斥力を表わす項であり，b∗ij · fA (rij ) は引力を表わす項である．またfc (rij ) はカットオフを表わす項である．bij は，原子 i，原子 j 以外に原子 k も含めた
3 個の原子によって定められる項であり，次式で表わされる．
bij = (1 + βin ζijn )−1/2n
ζij =

∑

[

(2.21)]

fc (rik )g(θijk )exp λ33 (rij − rik )k(̸=i,j)[g(θijk ) = 1 + c2 /d2 − c2 / d2 + (h − cos θijk )2(2.22)](2.23)

ここで θijk は結合 rij ，rik 間の角度を示している (図 2.1)．
実際の計算では，ある原子 i からカットオフ距離 (R + D) 内に存在する原子をリストアップし，カットオフ距離内の原子 j, k, l, … に対し，結合 i − j, i − k, i − l, …の二体間に働く力，及びそれぞれの結合から二つ，例えば結合 i − j と結合 i − k を選
びその三体間によって発生する力をそれぞれ求め，これらの力のベクトルを重ね合わせることにより原子 i に働く力を求めた．

k

θ ijk
i

j

Fig.2.1 Three molecules i, j, k and bendig angle θijk

2.4 Tersoﬀ 型ポテンシャル

11

次にポテンシャルパラメータについて説明する．Si の場合には，T2(26) および T3(27)と呼ばれる 2 種類のパラメータセットが提案されている．T2 は，表面物性を優先させてフィッティングしており，クラスターの形状を正確に表現することができるが(34) ，エネルギーが ab-initio の結果と異なることが報告されている(35) ．一方，T3 は，T2 では実験と一致しない弾性定数 C44 を修正し，結合角が 109.47◦ になるようにフィッティ
ングしているため，バルクの記述には適しているが，結合角に依存した原子挙動を示すことが報告されている(34) ．Ga–Ga，As–As の結合パラメータには Smith がバルク
およびクラスターに合わせ込んだパラメータ(36) を，Ga–As の結合パラメータには M.

Sayed と J. H. Jeﬀerson が凝集エネルギーやバルク係数等をフィッティングしたパラメータ(37) を用いた．ポテンシャルパラメータを以下に示す．

Table 2.3 Parameters for silicon, gallium and arsenic

A [eV]
B [eV]−1
λ1 [Å ]−1
λ2 [Å ]−1
λ3 [Å ]βncdh
R [Å ]
D [Å ]Si − Si(T2)
3.2647×103
9.5373×101
3.2394
1.32581.32583.3675×10−1
2.2956×101
4.83812.04170.0
3.0
0.2

Si − Si(T3)
1.8308×103
4.7118×102
2.47991.7322
1.73221.1000×10−6
7.8734×10−1
1.0039×105
1.6217×101−5.9825×10−1
2.85
0.15

Ga − Ga
9.9388×102
1.36123×102
2.508421.49082
1.490822.35862×10−1
3.472907.62977×10−2
1.97965×101
7.14592
3.5
0.1

As − As
1.57186×103
5.46431×102
2.384131.72873
1.72873
7.48809×10−3
6.08791×10−1
5.273137.51027×10−1
1.52924×10−1
3.5
0.1

Ga − As
2.54330×103
3.14460×102
2.828091.72301
1.72301
3.57192×10−1
6.317471.22630
7.90396×10−1−5.1849×10−1
3.5
0.1

2.4 Tersoﬀ 型ポテンシャル12原子に働く力は式 (2.2) で示すようにエネルギーの空間勾配によって表わされる．
[

]

∂
∂
r ji
∂ ∗
∂Eij
Fi = −=−(fc (rij )fR (rij )) − b∗ij
(fc (rij )fA (rij ))− fc (rij )fA (rij )
b
∂r i
∂rij
∂rij
rij
∂r i ij
[
]
∂Eij
∂
r ij
∂ ∗
∗ ∂
Fj = −=−(fc (rij )fR (rij )) − bij(fc (rij )fA (rij ))− fc (rij )fA (rij )
b
∂r j
∂rij
∂rij
rij
∂r j ij∂ ∗∂Eij= −fc (rij )fA (rij )b(2.24)
Fk = −∂r k∂r k ij
ここで次の関係を用いている．

∂rij
∂rij
r ji ∂rik
∂rik
r ki
=−
=
,
=−
=
∂r i
∂r j
rij
∂r i
∂r k
rik

(2.25)

r ij = r j − r i(2.26)

F i ，F j の第一項は結合 i − j 間に働くの力である．これらの力は結合 i − j に関して独立に計算しておく．残りの項については結合 i − j と他の結合 (i − k, i − l, …) との三体間を考慮して， b∗ij 及びその勾配を求めた上で計算する． b∗ij の勾配について以下
に示す．

∂b∗ij
1
=
∂r m
2

(

∂bij
∂bji
+
∂r m ∂r m)(2.27)1∂bij
1∂= − (1 + β n ζijn )− 2n −1(1 + β n ζijn )
∂r m
2n
∂r m
∂
∂ζij
(1 + β n ζijn ) = nβ n ζijn−1
∂r m
∂r m



(2.28)(2.29)



∂ζij
∂  ∑
=
fc (rik )g(θijk ) exp(λ33 (rij − rik )3 )
∂r m∂r m k(̸=i, j)=∑∂
fc (rik )g(θijk ) exp(λ33 (rij − rik )3 )∂r mk(̸=i, j)(m = i，j ，k)

ここで ζijk = fc (rik )g(θijk ) exp(λ33 (rij − rik )3 ) とおくと，∂ζijk∂ζijk ∂rik∂ζijk ∂(cos θ) ∂ζijk ∂rij
=
+
+
∂r i
∂rik ∂r i∂(cos θ) ∂r i
∂rij ∂r i∂ζijk ∂(cos θ) ∂ζijk r ji
∂ζijk r ki++=
∂rik rik∂(cos θ) ∂r i
∂rij rij(2.30)

2.4 Tersoﬀ 型ポテンシャル
(

13

)

∂fc (rik )r ki=− 3λ33 (rij − rik )2 fc (rik ) g(θ) exp(λ33 (rij − rik )3 )∂rikrik∂g(θ) ∂(cos θ)+fc (rik ) exp(λ33 (rij − rik )3 )∂(cos θ) ∂r i
r ji
+3λ33 (rij − rik )2 fc (rik )g(θ) exp(λ33 (rij − rik )3 )
rij
∂ζijk
∂g(θ) ∂(cos θ)= fc (rik ) exp(λ33 (rij − rik )3 )∂r j∂(cos θ) ∂r j
r ij
+3λ33 (rij − rik )2 fc (rik )g(θ) exp(λ33 (rij − rik )3 )
rij
(
)
∂ζijk
∂fc (rik )r ik=− 3λ33 (rij − rik )2 fc (rik ) g(θ) exp(λ33 (rij − rik )3 )
∂r k
∂rik
rik
∂g(θ) ∂(cos θ)+fc (rik ) exp(λ33 (rij − rik )3 )∂(cos θ) ∂r k

(

)

(

∂(cos θ)cos θ
1 r ji
cos θ
1
= −
+
+ −
+
∂r i
rij
rik rijrikrij()
cos θ r ij
1 r ik∂(cos θ)= −+∂r j
rij rijrij rik()∂(cos θ)cos θ r ik
1 r ij= −+∂r k
rik rikrik rij)(2.31)

r ki
rik(2.32)

また，∂g(θ)
2c2 (cos θ − h)= 2∂(cos θ)(d + (h − cos θ)2 )2である．

(2.33)

2.5 速度スケーリング法

2.5

14

速度スケーリング法

分子動力学法で温度制御する場合，もっとも簡単で直接的な方法として速度スケー
リング法がよく用いられる．熱統計力学より系の運動エネルギー K は次のように表さ
れる．

K=

N
1∑
3
mi (v i ·v i ) = N kB T
2 i=1
2

(2.34)

ここで，mi は原子 i の質量，v i は原子 i の速度，N は系の全原子数，kB はボルツマン
定数，T は系の温度である．式 (2.34) より，系の温度 T は原子速度を用いて，次のよ
うに求められる．
∑T =mi (v i ·v i )3N kB(2.35)

設定温度が TC ，式 (2.35) より求めたある時刻の温度が T のとき，速度スケーリング√

法では，各原子の速度 v i を

TC /T 倍し設定温度 TC に近づける．ベルレ法では，∆r i (t + ∆t) = r i (t + ∆t) − r i (t) = r i (t) − r i (t − ∆t) + (∆t)2

F i (t)m(2.36)√

を

TC /T ∆r i (t + ∆t) で置き換えることに相当する．平衡状態では，能勢の方法(38) など外部との熱のやりとりをする変数を考慮した拡張系の分子動力学法によって得られるカノニカルアンサンブルに一致することが示されている．

2.6 高速化手法

2.6

15

高速化手法

領域分割による高速化
式 (2.15) からわかるように，N 個の原子からなる系では，Etot の評価に N ×(N −1) 回の原子対の計算が必要となる．一方，実際の結晶中では近接原子による遮蔽 (screening)
効果により第二近接距離程度より離れた原子はほとんど作用を及ぼさないことが知られている．このため，分子動力学計算では相互作用打ち切り (カットオフ) 半径 rc を導入し (図 2.2)，その半径内の原子からの寄与のみを考慮する．
しかしながら，相互作用する原子対の検索に N × (N − 1) 回の試行を要するため，系が大きくなるにつれ計算負荷が飛躍的に増加する．これを避けるために rc よりひとまわり大きい半径 rfc (図 2.2) 内の原子をメモリーに記憶し，rfc 内での原子対の探索とすることによりオーダー N の計算に近づける方法 (粒子登録法(38) ) がこれまでよく用いられてきた．しかしながら，粒子登録法では rfc 半径より外の原子が rc 内に達すると力の評価が適切でなくなるので，一定のステップ毎に登録粒子の更新 (N × (N − 1)回の探査) を行わなければならない．このため，系がある程度の規模以上に大きくなると，粒子登録による高速化は登録更新の負荷により打ち消される．

rfc

rc

Fig.2.2 Schematic of bookkeeping method

2.6 高速化手法

16

領域分割法では，まず図 2.3 に模式的に示すようにシミュレートする系をカットオ
フ距離程度の格子状に分割する．ある原子に作用する力を評価する際には，その原子が属する領域（図 2.3 の着色部）と隣接領域内 (図 2.3 の斜線部) の原子からカットオフ
距離内の原子を探索する．原子が属する領域は，位置座標を領域ブロックの辺長 bx , by
で除した際の整数により判断できるので，領域分割そのものの計算負荷は小さい．領域分割法は，粒子登録法において登録更新の負荷が大きくなるような大規模な系の高速化に適している．

y
by
0

x
bx

Fig.2.3 Schematic of domain decomposition

第3章

Ni 基超合金 γ/γ ′ 界面における転位挙動
3.1

シミュレーション条件

3.1.1

解析モデル

Ni 相の原子面数に対して，一原子面少ない Ni3 Al 相を接合して半整合界面を作成した (図 3.1(a))．Ni3 Al 相は 49 × 49 × 35 の L12 格子を含み，Ni 相は 50 × 50 × 35 の fcc
格子を有する．原子数は 686140 個である．初期配置作成時にはミスフィットが周期境界に集中しないよう，Ni 相の格子定数 aγ を 0.3498 nm，Ni3 Al 相の格子定数 aγ ′ を0.3570 nm とし，50aγ = 49aγ ′ としている．このセルの x, y 方向に周期境界を適用し，z 方向は自由表面の境界条件で 15000 fs の初期緩和計算を行った．このとき，x, y 方向の応力が零となるようにセル寸法を制御している．また，熱揺動の影響を排除するため
温度は 10 K とし，速度スケーリングにより制御した．初期緩和後に得られたミスフィッ
ト転位を図 3.1(b) に示す．図では γ ′ 相の原子，および，Common Neighbor Analysis(CNA)(39) により fcc でも hcp でもないと判断された原子 (転位芯) のみ表示している．
また，後述のシミュレーションで導入する転位のすべり面との関係から，[110] 方向の転位線ベクトルを t⊥ ，[1̄10] 方向のそれを t// と表記することにする．それぞれの転位のバーガースベクトルは b⊥ = a/2[1̄10], b// = a/2[1̄1̄0] であり，基本的に転位線に垂直 (刃状転位) である(12) ．ここで a は平均格子長さである．

17

3.1 シミュレーション条件

50 fcc lattices(Ni)(17.6nm)γ

35 lattices(12.3nm)fcc

γ'

35 lattices(12.6nm)50

18

49 L1 la2 ttices(Ni3Al) (17.6nm)49

L1 2(a) Cell dimensionst//t⊥ b⊥ b//y [010]z [001]x [100](b) Line and Burgers vectorsof misfit dislocation

Fig.3.1 Simulation cell and misﬁt dislocations

3.1.2

刃状転位シミュレーション

図 3.2 に模式的に示すように，周期境界を考慮して Ni 相の上側表面を domain1 とdomain2 に分割し，どちらか一方を押し込むことで [1̄10] 方向に無限に長い刃状転位を導入するシミュレーションを行った．domain1 と domain2 の境界にある (111̄) 面は，
ミスフィット転位と平行または垂直に交わるすべり面である (図 3.3 左)．domain1 を押し込んだときと domain2 を押し込んだときでは，ミスフィット転位に平行・垂直に交わる刃状転位のバーガースベクトルが逆になる．押し込みは，図 3.3 右に示したように (111̄) 面上の 2 つの部分転位による原子移動 b1 , b2 に分けて，それぞれ 2500 fs で導入されるように上端面原子に微小変位を与えて制御した．導入した転位を界面に接近させるために，上端面に変位を与え続けて次々と転位を発生させた．

3.1 シミュレーション条件

Domain 1 Domain 2(a) Indentation 1
[001]
[100]
[010]

Domain 1

Domain 2x [100](b) Indentation 2

y [010]

Fig.3.2 Indentation for injection of edge dislocations

19

3.1 シミュレーション条件

20

Atoms on the upper slip plane
Atoms on the lower slip planeab1 = [11 2]
6ab = [1 0 1]2ab2 = [2 11]
6

Fig.3.3 Slip plane and displacement control for injection of edge dislocations

3.1.3

らせん転位シミュレーション

刃状転位シミュレーションと同様に，周期境界を考慮して図 3.4(a) のように上端面を 2 つの領域に分ける．domain2 は 2 つの部分転位による原子移動 b3 , b4 に分け，それぞれ 2500 fs かけて微小変位を与えた (図 3.4(b) 左下の原子移動)．しかしながら，単純に片方の domain を [1̄10] 方向にずらしても上端面の原子が移動するだけであり，らせん転位が結晶内部に進行することはない．そこで，図 3.4(b) 右に示すように，転位線に対して垂直方向の変位を段階的に domain1 に与えることで (111̄) 面上の [1̄10] 方向
らせん転位を導入した．domain1 のうち，図で左上側の境界に近い部分 (色の濃い部
分) では domain2 の b3 の変位を与える際に b5 = a/6[112] だけ移動させ，b4 の変位導入時には b6 = a/6[1̄1̄2̄] の移動を与えて元に戻す (図 3.4(b) 右上)．domain1 内でこの
垂直方向変位を段階的に少なくしていき，図中右下側の境界 (色の薄い部分) では零とした．以上の変位制御により，domain1 から見て右下側の境界では domain2 が b3 +b4
だけ相対変位するが，左上側の境界では domain2 から見ると domain1 が b7 +b8 だけ
相対変位する (図 3.4(c))．

3.1 シミュレーション条件

Domain 1

21

Domain 2
Domain 1

Domain 2

[100]
[001]
[100]

[010]

[010]

(a) Domains on upper surfaceab5 , b6 = ± [112]
6

b4

a
a
b7 = [2 11]
b4 = [2 11]6aa6b3 = [1 2 1]
b8 = [1 2 1]
66b3(b) Atom migration control(c) Resulting atom migration

Fig.3.4 Displacement control for injection of screw dislocations

3.1.4

せん断シミュレーション

図 3.5(a) に模式的に示すように，Ni 相の上端面原子を chuck1，Ni3 Al 相の下端面原子を chuck2 とし，それぞれに反対方向の変位を与えてせん断シミュレーションを行った．制御対象となる原子は chuck1 および chuck2 ともに厚さ 1 格子分であり，chuck1
は 0.350 nm，chuck2 は 0.357 nm である．本解析では [100] 方向および [110] 方向にせん断をかける 2 つのシミュレーションを行った．[100] 方向へのせん断シミュレーショ
ンでは，chuck1 に対して毎ステップ [100] 方向に dux だけの微小変位を与え，chuck2
に対して毎ステップ [1̄00] 方向に dux だけの微小変位を与えた．ここで dux は，

1
dux = zz {tan(γ + dγ) − tan(γ)}
2

(3.1)

3.1 シミュレーション条件

22

Chuck 1

dux

γ

zz

dγ

zz
[001]

Chuck 2

[001]
[010]

[100]

[100]

(a) Controlled domain(b) Shear displacement in [100]Fig.3.5 Shear simulation procedure

である (図 3.5(b))．dux は制御部分の変位量，γ はせん断ひずみである．ひずみ速度 dγを一定とし，dγ = 1.0 × 10−6 [ 1 / fs ] = 1.0 × 109 [ 1 / s ] とした．したがって dux は，式(3.1) に基づいて毎ステップ変化させている．一方，[110] 方向へのせん断シミュレー
ションでは，[100] 方向へのせん断ミュレーションと同じく，上下端面をそれぞれ制御
原子とし，chuck1 に対して毎ステップ [110] 方向に dux だけの微小変位を与え，chuck2
に対して毎ステップ [1̄1̄0] 方向に dux だけの微小変位を与えた．

3.2 シミュレーション結果および考察

3.2

シミュレーション結果および考察

3.2.1

刃状転位シミュレーション 1

23

domain1 を押し込んだ際の転位挙動を CNA により可視化し，界面近傍の状態を図 3.6に示した．t = 9000 fs の図 3.6(a) において示したように，上端面より発生した転位は薄く着色した leading partial と trailing partial の芯の間に，濃く着色した積層欠陥として存在している．ミスフィット転位と平行に衝突する転位を転位 I，ミスフィット転

DislocationⅡ

DislocationⅠbⅠbⅡ
Trailing partial
Leading partial(a) t = 9000 fs(b) t = 10000 fs(c) t = 11000 fs

b Ⅰ' b Ⅰb//(d) t = 12000 fs(e) t = 13000 fs(f) t = 14000 fs(g) t = 15000 fs(h) t = 16000 fs(i) t = 17000 fs

Fig.3.6 Snapshots of dislocations in indentation of domain 1

3.2 シミュレーション結果および考察

24

位と垂直に切り合う転位を転位 II とする．図 (a) 中に記号で示したように，押し込みによってこれらのすべり面の間に余剰原子面が導入される．転位 II が転位 I より広く拡張しているのは，自由表面から界面までの距離の差によるものと考えられる．紙面
手前方向にそれぞれの転位線ベクトルをとれば，転位 I，II のバーガースベクトルはbI = −(b1 +b2 ), bII = −bI となる．
まず，ミスフィット転位に平行に衝突する転位 I について説明する．この転位は t =12000 fs の図 (d) で界面に到達する．図中，右と左の界面転位にほぼ同時刻に達しているが，これらは周期境界によるイメージではなく同一の転位線の別の部分である．いずれも，leading partial がミスフィット転位に衝突し動きを止める．図 3.7 に 11700 ∼ 12900 fs
における衝突の様子を詳細に示すように，ミスフィット転位の位置で止められた 1 本目の転位は，後続の転位の力を受けて拡張転位の幅が小さくなるとともに，leading partialは (111̄) すべり面に交差する (111) すべり面上で，[1̄1̄2] 方向ならびに [112̄] 方向に広がるように移り変わる (図 3.7(c))．最終的に t = 12600 fs の図 3.7(d) では (111) 面上の拡
張転位となり，界面から遠ざかるように進行した (図 3.6(f)-(i))．結果的に転位 I は反射するように界面で跳ね返った．跳ね返った後の転位 I の転位線ベクトルを紙面手前方向に定義すれば，反射した転位 I のバーガースベクトル b′I は転位 I のバーガースベ
クトル bI とミスフィット転位のバーガースベクトル b// の和である (図 3.6(f) の矢印の3D view[11 2][11 2]Side view(a) 11700 fs(b) 12000 fs(c) 12300 fs(d) 12600 fs

Fig.3.7 Reﬂection of dislocation I(e) 12900 fs

3.2 シミュレーション結果および考察

25

関係)．その後，転位 I は界面から離れるが，図 3.6(h) および (i) から分かるように，界面上のミスフィット転位が消滅している．また，t = 14000 ∼ 15000 fs にかけて，2 本
目の転位 II ′ の trailing partial と左側の界面転位で反射した転位 I の leading partial が合体し，再び (111̄) 面にすべり系を変えているのも興味深い．
次に，転位 II がミスフィット転位と垂直に切り合う際の挙動について検討する．t =10000 fs の図 3.6(b) を見ると，転位 II の leading partial がミスフィット転位に引き寄せられて湾曲している．その後，後続の転位 II ′ の力を受けて転位の幅が小さくなるとともに転位の湾曲も小さくなる (図 (c), (d))．t = 13000 ∼ 17000 fs にかけて転位 II は

2 本のミスフィット転位によりせん断されるようにずれが生じた．t = 15000 fs の転位
を上から見たものを図 3.8 に示す．図中矢印で示したように，ミスフィット転位も相対的にずれている．この切り合いの様子を詳しく観察するため，t = 14000 fs における転位の構造を，図 3.9 に拡大して示す．図 (a) に示したように様々な角度から見た切り合いの様子を図 (b)〜(d) に示している．図 (c) から分かるように，ミスフィット転位から (1̄11) および (11̄1) すべり面上に leading partial が発生している．この部分転位は
図 3.9(a) に示すように視点 (c), (d) 側にのみ発生し，視点 (b) 側には見られない．転位

II はこのくさび状の拡張転位の外側を進行していることが図 (d) より分かる．図 3.10(a)にこれらのくさび状の拡張転位のすべり面と衝突させた刃状転位のすべり面の関係を，
図 3.10(b) に図 3.9(c) の転位群のバーガースベクトルの関係を模式的に示す．図 3.10(a)に示したように，(111̄) 面上の転位 II による原子移動は b1 +b2 = a/2[1̄01̄] である．一

Fig.3.8 t = 15000 fs，topview

3.2 シミュレーション結果および考察

26

方，(1̄11) 面上の leading partial のバーガースベクトルは a/6[1̄12̄] であるが，a/6[2̄1̄1̄]の trailing partial によって生じる原子移動は a/2[1̄01̄] となり (111̄) 面上の転位 II のそれと一致する． 転位 II の leading partial のバーガースベクトル b1 = a/6[1̄1̄2̄] と，ミスフィット転位から (1̄11) 面に生じた trailing partial のそれ (a/6[2̄1̄1̄]) は一致しないが，すべての成分が負であり方向が近い．転位 II の trailing partial のバーガースベクトル

b2 = a/6[2̄11̄] と (1̄11) 面の leading patial のそれ (a/6[1̄12̄]) の関係も同様である．このため，図 3.10(b) の左側ではこれらの転位線の交点が (111̄)，および (1̄11) 面の交線上でスライドすることにより bII の転位が γ ′ 相に侵入することが可能である．一方，(11̄1)
面側に生じる leading partial のバーガースベクトルは a/6[11̄2̄] で，転位 II の leading

partial (b1 = a/6[1̄1̄2̄]) とジャンクションを形成しながら拡張しているが，バーガース
ベクトルが a/6[1̄2̄1̄] の trailing partial を生じても拡張転位のそれは a/2[01̄1̄] となり転位 II による原子移動と一致しない．このため，これ以上拡大して γ ′ 相に侵入することはなかった．

[0 0 1](c)[11 0][11 0](d)(b)(a) View point

DislocationⅡ
Misfit dislocation(b) Rotated view 1
Misfit dislocation[0 0 1][0 0 1][11 0][11 0]

Misfit dislocation(c) Rotated view 2[11 0][11 0](d) [110] viewFig.3.9 Magniﬁed view of dislocations in Fig.4.1(f)3.2 シミュレーション結果および考察

abⅡ = [1 0 1]
2

(111)

a
b1 = [11 2]6ab 2 = [2 11]6ab ⊥ = [11 0]2
a[11 2]
6a[2 11]
6

(111)
bⅡ

(111)

[001]
[100]
[010]

(a) Slip planesab ⊥ = [11 0]2ab 2 = [2 11]6a[2 11]
6abⅡ = [1 0 1]
2

b2

bⅡ

a
b1 = [11 2]
6

b⊥
bⅡ

a
b1 = [11 2]
6a[11 2]
6a[11 2]
6(b) Burgers Vectors
Fig.3.10 Schematic of slip planes and Burgers vectors in Fig.4.4(c)273.2 シミュレーション結果および考察

3.2.2

28

刃状転位シミュレーション 2

domain2 の押し込み下で生じる転位挙動を図 3.11 に示す．前項と同様，ミスフィット
転位に平行に衝突するものを転位 I，垂直に交差するものを転位 II と称することにする．
まず，転位 I の挙動について説明する．転位 I の leading partial がミスフィット転位に到達した t = 12000 fs の図 (b) では，ミスフィット転位から γ ′ 相側に leading partial が発生している．その leading partial は γ ′ 相中に大きく張り出すと共に，交差するミス
フィット転位から，前項で観察したような (1̄11) および (11̄1) 面への leading partial の発生が認められる (図 (c)-(g) 左側のミスフィット転位部分)．その後，ミスフィット転位の交差部分から穴が開くように転位が分離し，その穴が拡大するようにして前縁の拡張転位が γ ′ 相内に侵入した (図 (g)-(j) 丸で囲った部分)．図 3.12 に t = 15000 ∼ 18000 fs
におけるその挙動を拡大して示す．図 (a) は観察した転位の位置と視点の向き，図 (b)〜(e) がその時間変化である．t = 15000 fs の図 (b) において，上方の転位はミスフィッ
ト転位に接近している後続の転位 I ′ である．下方の幅の広い転位は，ミスフィット転位から発生した leading partial，およびミスフィット転位と衝突した転位 I の trailing

partial からなる．また，前項で述べた交差すべり面への leading partial の発生が認められる (矢印 ,)．ミスフィット転位の交差部分は多数の転位の交点となり，fcc でもhcp でもないと判別された原子が多くなっている．t = 16000 fs の図 (c) において，欠陥原子が集中した交差部分に小さな穴が発生し，それが拡大していく様子が図 (d), (e)から分かる．また，図 (d), (e) で一番下側の転位は，ミスフィット転位から発生した(1̄11) 面上の leading partial との交差部分をスライドさせながら γ ′ 相へ侵入しており(図中丸で囲った部分)，図 3.10(b) で述べたメカニズムと同じである．一方，界面上にはミスフィット転位が残されているが，最初と同じ位置にはなく界面上でずれている．
その様子を詳しく観察するため，界面上下 1.25 nm の薄板領域の転位挙動を上から見たものを図 3.13 に示す．図では周期境界を含めた全体的な形態変化を示すために，セルを 2 つ並べて示している．図の中央部分が先述した転位の交点であるが，図 (c) に矢印で示したように同じバーガースベクトル・すべり面の関係にある別の交差部にもループ状の転位が確認できる．

3.2 シミュレーション結果および考察

bⅡbⅠ(a) t = 11000 fs(b) t = 12000 fs(c) t = 13000 fs(d) t = 14000 fs(e) t = 15000 fs(f) t = 16000 fs(g) t = 17000 fs(h) t = 18000 fs(i) t = 19000 fs(j) t = 20000 fs(k) t = 21000 fs(l) t = 22000 fs

Fig.3.11 Snapshots of dislocations in indentation of domain 2

29

3.2 シミュレーション結果および考察(a) View point①②(c) t = 16000 fs(b) t = 15000 fs(d) t = 17000 fs(e) t = 18000 fs

Fig.3.12 Dissociation from the junction between misﬁt dislocation node andapproaching edge dislocation

Ⅱ
Ⅰ

Ⅰ
Ⅱ
Ⅰ

(a) t = 16000 fs(b) t = 17000 fs(d) t = 19000 fs(e) t = 20000 fs(c) t = 18000 fs

Fig.3.13 Morphology of dislocations on the interface (Indentation 2)30

3.2 シミュレーション結果および考察

31

次に転位 II の挙動について説明する．t = 13000 fs の図 3.11(c) において垂直なミス
フィット転位と接触するが，刃状転位シミュレーション 1 のようにミスフィット転位から引き寄せられることはなく大きな湾曲も見られない．その後，後続の転位 II ′ の力を受けて γ ′ 相に侵入しようとするが，t = 20000 fs の図 (j) から分かるように，転位 IIはミスフィット転位と切り合って侵入することはなく，ミスフィット転位にピンニング
されている．刃状転位シミュレーション 1 との違いを明確にするために，[110] 方向に見た転位 II とミスフィット転位の交差部を拡大して図 3.14 に示した．図より，刃状転位シミュレーション 1 では転位 II の交差部分が収縮しているのに対して，刃状転位シ
ミュレーション 2 における転位 II は拡張したまま堆積している．また，交差するミス
フィット転位から別のすべり面への leading partial の発生も認められない．
刃状転位がミスフィット転位と切り合って γ ′ 相へ侵入したのはシミュレーション 1
における転位 II，シミュレーション 2 における転位 I であり，いずれも直交するミス
フィット転位から交差するすべり面へ leading partial を生じていた．どちらも，すべり面に対して上側に余剰原子面があり，ミスフィット転位の余剰原子面を考えると，傾いてはいるが同符号の転位と考えられる (図 3.15(a))．一方，すべり面に対して下側に余剰原子面がある場合 (図 3.15(b))，刃状転位シミュレーション 1 の転位 I は，異符号で平行なミスフィット転位と合体して別の転位となった (図 3.6(f) の矢印の反応) が，シミュレーション 2 の転位 II では垂直に切り合うため，交差部が不動転位となりピンニングされたものと思われる．

3.2 シミュレーション結果および考察(a) Indentation 1 (b) Indentation 2
Fig.3.14 Magniﬁed view of edge dislocations crossing over misﬁt dislocations

Ⅱ

Ⅰ

Cross the line

On the line(a) Same sign
Ⅰ

Ⅱ

On the line

Cross the line(b) Opposite sign

Fig.3.15 Schematic explanation for interaction between misﬁt dislocation andedge dislocation32

3.2 シミュレーション結果および考察

3.2.3

33

らせん転位シミュレーション

上端面よりらせん転位を導入したシミュレーションにおける転位の挙動を図 3.16 に示す．これまでと同様，ミスフィット転位に平行に衝突するすべり面の転位を転位 I，
垂直に交差するものを転位 II と称することにする．まず，ミスフィット転位に平行に衝突する転位 I について説明する．t = 12000 fs の図 (c) において転位 I がミスフィッ
ト転位に接触すると，ほぼ同時に交差すべり面への leading partial が発生し，鏡のよ
うに跳ね返る．刃状転位シミュレーション 1 では，後続の転位による拡張転位の圧縮
後，交差すべり面への乗り換えを生じていたのに対し，らせん転位では極めてスムー
ズに跳ね返っている．図 3.17 に，刃状転位シミュレーション 1 におけるすべり面の乗り換えと，らせん転位シミュレーションにおける挙動の違いを横から見て比較している．らせん転位の場合，γ ′ 相に入ることなく反射するようにすべり面を変えていることが分かる．一方，跳ね返った転位 I の leading partial は (111) 面上を [1̄1̄2] 方向に進行するが，trailing partial が界面から離れる前に 2 本目の転位 II の leading partial と合体する (図 3.16(g) 矢印)．跳ね返った転位 I のバーガースベクトルと 2 本目の転位 II のそ
れは符号が逆であるため対消滅し，拡大した積層欠陥を狭めるように転位 I の leading

partial は戻って行く (図 (h), (i))．
次に，ミスフィット転位と転位 II の切り合いについて述べる．t = 17000 fs の図 (h)で界面に到達した後，拡張転位の幅が小さくなるとともにミスフィット転位との切り合いを生じ，界面上のミスフィット転位の形状が変化する (図 (i)-(o))．この形状変化
を図 3.18 に模式的に示す．切り合いによって

ミスフィット転位，

らせん転位，

隣のミスフィット転位，とバーガースベクトルの向きが等しい転位群がジグザグに繋がって 1 本の転位を形成する (図 3.18(b))．この転位は直線状になるように界面上で変形する (図 3.18(c))．t = 21000 fs の図 3.16(l) および図 3.18(c) では，初期の転位網の対角線方向にミスフィット転位が配向している．さらに，らせん転位に平行なミスフィッ
ト転位上で転位の端点が移動し，最終的に転位網の位置が初期のそれより 1/2 だけずれる (図 3.16(o) および模式図 3.18(d))．さらに転位を導入し続けると，同様の切り合い，および，界面上の転位網の形態変形を生じ，転位網が 1/2 ずれて t = 33000 fs では
図 3.19 に示すように元に戻った．

3.2 シミュレーション結果および考察

bⅡbⅠ(a) t = 10000 fs(b) t = 11000 fs(c) t = 12000 fs(d) t = 13000 fs(e) t = 14000 fs(f) t = 15000 fs(g) t = 16000 fs(h) t = 17000 fs(i) t = 18000 fs(j) t = 19000 fs(k) t = 20000 fs(l) t = 21000 fs(m) t = 22000 fs(n) t = 23000 fs(o) t = 24000 fs

Fig.3.16 Snapshots of dislocation in shear simulation for screw dislocation

34

3.2 シミュレーション結果および考察(a) Indentation1(b) Shearing

Fig.3.17 Reﬂection of edge and screw dislocation at parallel misﬁt dislocation

③
①

②

(a)

(b)

(c)

(d)

Fig.3.18 Schematic of morphology change in misﬁt dislocation (shear simulation for screw dislocation)

Fig.3.19 Misﬁt dislocation network at t = 33000 fs (shear simulation for screw
dislocation)35

3.2 シミュレーション結果および考察

36

せん断シミュレーション

3.2.4[100] 方向にせん断をかけたときの転位挙動を図 3.20 に，せん断応力 τzx とせん断ひ
ずみ γzx の関係を図 3.21 に示した．図 3.20(a) は緩和後のミスフィット転位及び上下端面の制御部分である．t = 15000 fs の図 3.20(b) において，ミスフィット転位は界面上を [100] 方向にわずかに移動している．その後，ミスフィット転位はゆっくりとせん断
方向に動き出し，t = 31000 fs の図 3.20(c) において，ミスフィット転位はセルの半分だ
け [100] 方向に移動している．t = 36000 fs (γzx = 0.036) において，図 3.21 に示したよ
うにせん断応力は最大値 τzx = 3.08 GPa を示すが，その後は図 3.20(e)〜(g) に示した
ように，ミスフィット転位は点欠陥を残しながらほぼ一定速度で動き続けた．界面上の点欠陥は転位の通過に伴い消滅と発生を繰り返している．図 3.22 に t = 90000 fs におけるシミュレーションセル全体のスナップショットを示す．γ 相と γ ′ 相内には乱れが見られず (図 3.20(h))，界面の相対的なすべりで変形を吸収している (図 3.22)．

[0 0 1][1 1 0][1 1 0](a) t = 0 fs, γ zx = 0(b) t = 15000 fs, γ zx = 0.015(c) t = 31000 fs, γ zx = 0.031(d) t = 36000 fs, γ zx = 0.036(e) t = 40000 fs, γ zx = 0.04(f) t = 45000 fs, γ zx = 0.045(g) t = 50000 fs, γ zx = 0.05(h) t = 90000 fs, γ zx = 0.09

Fig.3.20 Motion of misﬁt dislocation under the [100] shear on the interface

3.2 シミュレーション結果および考察

Shear stress, τzx [GPa]43.08 GPa
3

2

1

0

0.02

0.04

0.06

0.08

Strain, γzx
Fig.3.21 Shear stress (τzx ) against strain under the [100] shear of the interface[0 0 1][0 1 0][0 1 0][ 1 0 0](a) Side view

[0 0 1][ 1 0 0](b) Rotated viewFig.3.22 Analytical cell at t = 90000 fs, γzx = 0.09

37

3.2 シミュレーション結果および考察

38

次に [110] 方向にせん断をかけたときの転位挙動を図 3.23 に示す．t = 15000 fs の
図 3.23(b) において，[100] 方向へのせん断シミュレーションと同じく，ミスフィット
転位は界面上を移動しはじめる．図 (a) において黒い矢印で示したミスフィット転位(せん断方向に垂直) が，図 (b) において図の右上方向 ([110] 方向) に移動しているのが分かる．x，y 方向に周期境界を適用しているため，黒い矢印で示したミスフィット転位は，図の右上からセル外に出ると図の左下から現れる．[110] 方向へのせん断シミュ
レーションでは，白い矢印で示したミスフィット転位 (せん断方向に平行) は動かず，黒い矢印で示した部分のみがせん断方向に移動した．また，[100] 方向へのせん断シミュ
レーションでは，転位網の交点の後に多く点欠陥が発生したが，[110] 方向へのせん断
シミュレーションではほとんど点欠陥が観察されていない．そこで，本シミュレーショ
ンでは黒い矢印で示した転位の位置を容易に決定することができる．図 3.24 に転位の位置の変化を示す．せん断ひずみ γz′ x′ 約 0.015 までは転位はほとんど動いていないが，その後 γz′ x′ = 0.04 以降は一定速度となった．せん断応力 τz′ x′ とせん断ひずみ γz′ x′ の関係を図 3.25 に示す．ここで x′ y ′ z ′ 座標は x′ = [110]，y ′ = [1̄10]，z ′ = [001] である．
ミスフィット転位の移動速度が一定となる γz′ x′ = 0.04 においてせん断応力の増加が鈍化し，応力–ひずみ曲線の折れ曲がり点となっている．せん断応力はピークを持たず，ひずみ硬化を伴う弾塑性的な挙動を示している．先の [100] 方向へのせん断シミュレー
ションにおいて示された最大せん断応力は，ミスフィット転位が界面を滑るための「最大静摩擦力」であり，その後界面には多数の点欠陥など大きな構造緩和を生じるため応力が急減する．一方，[110] 方向へのせん断では，1 本の転位線 (周期的に並んでいるが) が一定速度で運動する動摩擦力に対応するものと考えられる．その後のせん断応力のゆるやかな増加は図 3.23(f)〜(h) において転位線が後方に残像のように乱れていることから，動点な効果と推測される．なお，図 3.20 及び図 3.23 において最大せん断応
力が異なるのは，結晶の異方性によるものである．

3.2 シミュレーション結果および考察[0 0 1][1 1 0][1 1 0](a) t = 0 fs, γ z'x' = 0(b) t = 15000 fs, γ z'x' = 0.015 (c) t = 25000 fs, γ z'x' = 0.025 (d) t = 35000 fs, γ z'x' = 0.035(e) t = 42500 fs, γ z'x' = 0.0425 (f) t = 50000 fs, γ z'x' = 0.05(g) t = 55000 fs, γ z'x' = 0.055 (h) t = 90000 fs, γ z'x' = 0.09

Fig.3.23 Motion of misﬁt dislocation under the [100] shear on the interface

Position, d, nm

80

60

40

16.1 nm

20

0

0.02

0.04

0.06

0.08

Strain, γz'x'
Fig.3.24 Position of the dislocation in Fig.3.23(a) against strain39

3.2 シミュレーション結果および考察

Shear stress, τz'x'[GPa]43.28 GPa
3

2

1

0

0.02

0.04

0.06

0.08

Strain, γz'x'
Fig.3.25 Shear stress (τz′ x′ ) against strain under the [100] shear of the interface40

3.3 結言

3.3

41

結言

半整合界面上のミスフィット転位と単相中を進行する転位の転位芯同士の相互作用を明らかにするため，格子数の異なる γ, γ ′ 相を接合して作成した半整合界面上のミ
スフィット転位に，直線状の刃状またはらせん転位を平行または垂直に衝突させるシ
ミュレーションを分子動力学法により行った．また，外力下でのミスフィット転位の形態変化を明らかにするため，同モデルを用いて [100] 方向または [110] 方向へのせん
断シミュレーションを行った．得られた結果を要約して以下に示す．

I. ミスフィット転位と γ 相を進行する転位の衝突(1) 同符号側に余剰原子面をもつ刃状転位がミスフィット転位に衝突する場合，交差するミスフィット転位から交差すべり面へ等価なバーガースベクトルの転位を発生し，γ ′ 相内に侵入した．

(2) 異符号側に余剰原子面をもつ刃状転位はミスフィット転位と平行に衝突すると，合成して別の転位となり交差すべり面上を跳ね返った．一方，垂直に衝突した
場合は切り合うことができず，刃状転位はミスフィット転位によりピンニングされた．

(3) らせん転位をミスフィット転位に平行に衝突させた場合，やはり交差すべり面上を跳ね返ったが，刃状転位の場合は拡張転位が収縮した後にすべり面を変えたのに対して，らせん転位の場合は leading partial が界面で反射されるようにすべり
面を変えていた．

(4) らせん転位をミスフィット転位に垂直に衝突させた場合，交差部分で転位の切り合いを生じ， ミスフィット転位， らせん転位， 隣のミスフィット転位，とジグザグに繋がった 1 本の転位を形成した．その後，転位が直線状になるように形を変え，初期の転位網の正方形対角線方向に配向した．さらに交差するミス
フィット転位上で端点が移動し，転位網の位相が 1/2 だけずれた．

3.3 結言

42

II. せん断力によるミスフィット転位の運動(5) [100] 方向にせん断をかけた場合も [110] 方向にせん断をかけた場合も，ミスフィッ
ト転位はせん断方向にゆっくりと移動し，その後ある程度ひずみが大きくなると，ミスフィット転位は一定速度で動き続けた．

(6) [100] 方向へのせん断では，ミスフィット転位の通過後に界面上に多くの点欠陥が発生し，大きな構造緩和を生じて応力が急減した．[110] 方向へのせん断では，せん断応力はピークを持たず，ひずみ硬化を伴う弾塑性的な挙動を示した．
以上総合すると，異符号の刃状転位，および，らせん転位とミスフィット転位の切り合いは，バーガースベクトルの合成で説明できる．一方，同符号の刃状転位が衝突する際には 2 つのすべり面の交線にジャンクションを有するような転位を生成することが明らかになった．また，外力下でもミスフィット転位が (001) 界面をすべり面として運動することも示された．本結果を離散転位動力学のような上位スケールシミュレー
ションへ反映することで，界面転位網のカッティング抵抗や形態変化などの界面転位
網のメカニズム解明が期待される．

第4章

Si/GaAs 界面における転位挙動
本章では Si/GaAs 界面の安定構造ならびに欠陥生成を MD シミュレーションにより検討することを目的とする．Si 層には Tersoﬀ が提案したパラメータセット(26)–(28) を用い，GaAs 層には M. Sayed と J. H. Jeﬀerson が提案したパラメータセット(37) を用いるが，4.1.2 項において述べるように，界面での結合 (Si–Ga，Si–As 等) の取り扱いが問題となる．本章では，まず Si/GaAs 界面エネルギーを第一原理計算と比較した後，Si 基板上に生成した GaAs ハットクラスターおよび界面の構造安定性を，各パラメー
タセットで検証する．その後，Si/GaAs 界面の引張シミュレーションを行う．

4.1

界面エネルギーの評価

4.1.1

三軸等方膨張・収縮を受ける Si・GaAs 単相のエネルギー評価

Si/GaAs 界面のポテンシャルパラメータの取り扱いを，まず界面エネルギーの点から議論する．第一原理計算で Si/GaAs 界面エネルギーを評価するが，現時点では計算能力の制約から Si，GaAs の単相セルを積み重ね，界面平行方向には周期性を仮定した解析となる．緒論で述べたように，Si と GaAs の格子ミスフィットは 4%にもなるので，Si，GaAs とも無ひずみとすることはできず，Si 相，GaAs 相ともひずんだ状態となる．したがって，界面エネルギーを算出するためには，Si 単相および GaAs 単相のエネルギーのひずみに対する変化を求めておく必要がある．そこで，Si，GaAs それぞれ単結晶について，図 4.1 に示すように単位格子を用い，格子長さ la の変化に対するエ43

4.1 界面エネルギーの評価

44

ネルギー変化を，Kresse らにより開発された平面波基底ウルトラソフト擬ポテンシャ
ル法(40) に基づく第一原理バンド計算コード VASP(41) (Vienna Ab-initio Simulation

Package)，ならびに Si，GaAs それぞれの Tersoﬀ 型ポテンシャルで評価した．第一原
理計算における計算条件を表 4.1 に示す．

Si 単結晶 (図 4.1(a)) および GaAs 単結晶 (図 4.1(b)) において格子定数 la(または原子間距離) を変化させたときの系のエネルギー変化を図 4.2 に示す．図 (a) において示した T2 および T3 のグラフは，表 2.3 第 2 列目および第 3 列目に示した Si のパラメー
タセットを用いている．図 (b) に示した GaAs のグラフは，表 2.3 第 6 列目に示した

Ga–As のパラメータセットを用いている．図 4.2 においてグラフ中に示した ◦ 印は，エネルギーの安定点 (極小点) である．第一原理計算でも同じ位置で極小値をとっており，
格子長さが正しく再現されている．エネルギーの絶対値は Si，GaAs いずれも第一原理より高いが，MD シミュレーションにおいて重要となるのはエネルギーの絶対値ではなく空間勾配である．そこで図 4.2 に示す第一原理計算の結果をグラフ上側にシフトして，極小点を一致させたグラフを図 4.3 に示す．図 (a) に示すように，Tersoﬀ 型
ポテンシャルの値 (T2 および T3) と第一原理計算の値 (VASP) との差を dE ，極小点におけるエネルギーを基準とした第一原理計算のエネルギーを EVASP とする．このとき，dE/EVASP が ±10%未満の領域を図の下側に矢印で示す．Si および GaAs ともに極小点近傍の勾配 (曲線形状) はほぼ等しい．これは体積弾性率をフィットさせている以
上，当然の結果である．図 (a) の引張側で T2 と T3 で誤差が 10%未満の限界が異なるのはカットオフ距離 (T2，T3 の曲線が大きく変化する点) によるところが大きい．

4.1 界面エネルギーの評価

: Ga

: Si

: As

la
z

y

la
z

x

la
la

y

x

la
la

(a) Si model(b) GaAs modelFig.4.1 Si and GaAs unit lattices

Table 4.1 DFT calculation condition

Number of ions
Number of electrons
Number of bands
Cutoﬀ energy [eV]
Number of k points

Si
GaAs
Si/GaAs
8
8
16
32
32
64
32
32
64
188.27 180.38
188.273
3
8
8
10 × 10 × 5

45

4.1 界面エネルギーの評価

VASP
T2
T3

0

-20

VASP
GaAs

60
40

Energy [eV]

Energy [eV]

20

20
0
-20

-40
-40
4

5

6

7

4

Lattice parameter [Å]

5

6

7

Lattice parameter [Å](a) Si(b) GaAs

Fig.4.2 Relationship between potential energy and lattice parameter on Si
and GaAs single crystal

VASP
T2
T3

0

la = 4.64
la = 6.44-20

la = 6.14

dE
EVASP

VASP
GaAs

60
40

Energy [eV]

Energy [eV]2020la = 4.94
0

la = 6.74

-20

-40
-40
4

5

6

Lattice parameter [Å](a) Si

7

4

5

6

Lattice parameter [Å](b) GaAs

Fig.4.3 Comparison of energy curve of Fig.4.2

746

4.1 界面エネルギーの評価

4.1.2

47

界面におけるポテンシャルパラメータの扱い

Tersoﬀ 型ポテンシャル(26)–(28) では，2 種類以上の元素を用いたシミュレーションを行う場合，式 (2.17)，(2.18)，(2.20) におけるパラメータ A，B ，λ1 ，λ2 ，R，D (以下：A〜D) を以下のように各元素のパラメータを相加・相乗平均して用いる．
Aij = (Ai Aj )0.5 , Bij = (Bi Bj )0.5 , Rij = (Ri Rj )0.5 , Dij = (Di Dj )0.5
λ1 ij = 0.5 (λ1 i + λ1 j ) , λ2 ij = 0.5 (λ2 i + λ2 j )(4.1)

一方，結合角に関する式 (2.21)〜(2.23) のパラメータ λ3 ，β ，n，c，d，h は，結合角を構成する周囲の原子には依存せず，中心の原子 i にのみに依存するとしている．例
えば，Si 原子と C 原子が結合している場合において，その 2 原子に働く結合エネル
ギーを求める際，パラメータ A は Si のパラメータ ASi と C のパラメータ AC を用いて

A = (ASi AC )0.5 となる．一方，M. Sayed と J. H. Jeﬀerson(37) が GaAs について提案した Tersoﬀ 型ポテンシャルでは，パラメータ λ3 ，β ，n，c，d，h を Ga–As 間について
フィッティングしており (表 2.3 第 6 列目)，中心原子 i だけでなく周囲の原子 j の種類にも陽に依存する．したがって，Si/GaAs 界面を扱う際に，これらのパラメータの扱いが問題となる．
パラメータ β ，n に関しては，式 (2.21) に示すように 2 体間の項に含まれるので，原子 i および原子 j の組み合わせが Ga 原子と As 原子の組み合わせであれば，Jeﬀersonの Ga–As のパラメータを用いるが，それ以外 (Si–Ga，Si–As 等) は原子 i に依存させる．一方，パラメータ λ3 ，c，d，h (以下：λ3 〜h) に関しては，式 (2.22)，(2.23) に示
すように 3 体間の項に含まれるので，以下に示す 2 通りの場合を考慮する．

[ 1 ] 原子 i，j および i，k の組み合わせが Ga 原子と As 原子の組み合わせであれば，
Jeﬀerson の Ga–As のパラメータを用い，それ以外は原子 i に依存させる[ 2 ] 原子 i，j の組み合わせが Ga 原子と As 原子の組み合わせであれば，Jeﬀerson のGa–As のパラメータを用い，それ以外は原子 i に依存させる
例えば，図 4.4 に示すような原子の組み合わせを考える．[ 1 ] のような取り扱いをする場合，図 (a) 及び (b) については原子 i と j の組み合わせ及び原子 i と k の組み合

4.1 界面エネルギーの評価

48

わせが Ga 原子と As 原子の組み合わせであるので Jeﬀerson の Ga–As の組に対するパラメータを用いるが，図 (c) のように Si 原子が混ざった場合や図 (d) のように同原子での結合が見られた場合は λ3 〜h を原子 i だけに依存するパラメータを用いる (すなわち，図の例では (c) は Ga–As の λ3 〜h ではなく，Ga 単体の λ3 〜h を用い，(d) では As
単体の λ3 〜h を用いる)．一方，[ 2 ] のような取り扱いをする場合，図 (a)〜(c) では原子 i と j の組み合わせが Ga 原子と As 原子の組み合わせであるので Jeﬀerson の Ga–Asのパラメータを用いるが，図 (d) の場合のみ，原子 i と j の組み合わせが原子 Ga と原子 As の組み合わせでないので，λ3 〜h は原子 i のパラメータ，つまり図の例では As のパラメータ (表 2.3 第 5 列目) を用いる．
なお，パラメータ A〜D に関しては，Tersoﬀ も Jeﬀerson も原子 i，j に依存させるとしているため，Si–Ga または Si–As の結合を考える場合には，式 (4.1) に従って相加・相乗平均により求めたパラメータを，Ga–As の結合を考える場合には，Jeﬀerson が提案したパラメータ (表 2.3 第 6 列目) を用いる．
以上のように，λ3 〜h に関して 2 通りのパラメータの取り扱い方がある．また，Siのパラメータセットには T2 および T3 の 2 通りがあり，この組み合わせも考えると合計 4 通りの取り扱い方がある．したがって，次項以降のシミュレーションにおいて以下に示す 4 通りの組み合わせについてそれぞれ検討する．

(1) Si のパラメータには T2 を用い，λ3 〜h に関しては [ 1 ] のような取り扱いをする(2) Si のパラメータには T2 を用い，λ3 〜h に関しては [ 2 ] のような取り扱いをする(3) Si のパラメータには T3 を用い，λ3 〜h に関しては [ 1 ] のような取り扱いをする(4) Si のパラメータには T3 を用い，λ3 〜h に関しては [ 2 ] のような取り扱いをする

4.1 界面エネルギーの評価

Ga i

As

Ga

Si

Ga

k

k

k

k

j As

(a) As-Ga-As

As i

j Ga

(b) Ga-As-Ga

Ga i

j

As

(c) As-Ga-Si

Fig.4.4 Combination pattern of atoms

As i

j

As

(d) As-As-Ga

49

4.1 界面エネルギーの評価

4.1.3

50

Tersoﬀ 型ポテンシャルによる界面エネルギーの評価

解析には，Si と GaAs を 1 格子ずつ接合した図 4.5 のスーパーセルを用いた．原子数は 16 個である．第一原理計算における計算条件は表 4.1 に示した．ここでは，Si 部，

GaAs 部の格子長さの差を考慮せず，同じ格子長さ la で格子点位置に Si，GaAs 原子を配置し，第一原理計算ソフト VASP ならびに 4.1.2 項の (1)〜(4) の組み合わせによる Tersoﬀ 型ポテンシャルで，セル辺長 la を変化させてエネルギーを評価した．

la を等方的に膨張収縮させたときの系のエネルギー変化を図 4.6 に，第一原理計算の結果をシフトして，極小点を一致させたものを図 4.7 に示す．図 4.6，図 4.7 共に，図(a) には界面での λ3 〜h の扱いを [ 1 ] としたポテンシャル (1) 及び (3) について，図 (b)には [ 2 ] の取り扱いのポテンシャル (2) 及び (4) について得られた結果を示している．
また，図 4.7 には，4.1.1 項と同じく dE/EVASP が ±10%未満となる領域を図下側に矢印で示している．ポテンシャル (1)〜(4) のいずれにおいてもエネルギーの安定点 (極小点) ならびにその近傍でのエネルギー勾配はほぼ一致している．

2la

z

y

x

la
la

Fig.4.5 Si/GaAs model4.1 界面エネルギーの評価

100

100

VASP
Potential (1)Potential (3)50Energy [eV]

Energy [eV]

50

0

-50

-100

VASP
Potential (2)Potential (4)

0

-50

4

5

6

-100

7

4

5

6

Lattice parameter [Å]Lattice parameter [Å](a) Potential (1) and (3)(b) Potential (2) and (4)7

Fig.4.6 Relationship between potential energy and lattice parameter on
Si/GaAs model

100

50

VASP
Potential (1)Potential (3)0la = 6.74

la = 4.94

la = 6.44
-50

-100

VASP
Potential (2)Potential (4)50Energy [eV]

Energy [eV]100la = 4.64

la = 6.44

0

la = 6.14
-50

4

5

6

7

-100

4

5

6

Lattice parameter [Å]Lattice parameter [Å](a) Potential (1) and (3)(b) Potential (2) and (4)Fig.4.7 Comparison of energy curve of Fig.4.6

751

4.1 界面エネルギーの評価

52

次に，格子定数の変化に対する界面エネルギーの変化を算出する．Si/GaAs モデルでは，Si 同士の結合，GaAs 同士の結合が存在し，系のエネルギーにはそれらのひずみエ
ネルギーも含まれる．そこで，界面だけのエネルギー変化を検討するために，Si/GaAs
モデルから 4.1.1 項で求めた Si 単相，GaAs 単相のエネルギーを引くことで界面エネル
ギーを算出する．Si 単相，GaAs 単相，Si/GaAs モデルの系のエネルギーをそれぞれE(Si)，E(GaAs)，E(Si/GaAs) とすると，界面エネルギー Eint は，Eint = [E(Si/GaAs) − E(Si) − E(GaAs)] /(2A)(4.2)となる．ここで A はスーパーセルの断面積である．
式 (4.2) により得られた界面エネルギーの，スーパーセルの la に対する変化を図 4.8に示す．図に示した範囲は，図 4.7 において示した dE/EVASP が ±10%未満の領域よりも 0.3 Å だけ大きい範囲としている．図中の ◦ 印は図 4.6 の各ポテンシャルの扱いならびに第一原理計算においてエネルギー極小となる la である．界面における λ3 〜h の扱
いを [ 1 ] とすると，格子定数の変化に対する界面エネルギーの変化が極めて急になり，第一原理計算の傾向と大きく異なる (図 (a))．一方，図 4.8(b) の破線のポテンシャル(2) は，la に対する変化は第一原理のそれに似ているが，界面エネルギーは負となっている (図 (b))．つまり，バルクよりも界面が存在する方が安定となり現実的でない．し
たがって，界面エネルギーの観点からは，ポテンシャル (4) の取り扱いが良いと思われる．すなわち，Si のパラメータには T3 を用い，界面で Si，Ga，As が乱れた部分では λ3 〜h の取り扱いを 4.1.2 項の [ 2 ] とする．

4.1 界面エネルギーの評価

VASP
Potential (1)Potential (3)

0.3

0.2

0.1

0

-0.1

5

6

VASP
Potential (2)Potential (4)0.4Interfacial energy [eV/ Å2]

Interfacial energy [eV/ Å2]

0.4

7

0.3

0.2

0.1

0

-0.1

5

6

Lattice parameter [Å]Lattice parameter [ Å](a) Potential (1) and (3)(b) Potential (2) and (4)7

Fig.4.8 Relationship between interfacial energy and lattice parameter on
Si/GaAs model53

4.2 界面原子構造の評価

4.2

54

界面原子構造の評価

前節では，界面エネルギーを評価するために Tersoﬀ 型ポテンシャルと第一原理計算を比較したが，扱った解析モデルが小さいため，Si 部，GaAs 部の格子長さの差を考慮せずに整合させた Si/GaAs 界面についてのものであった．本章では，第一原理計算では扱えない大きさの解析モデルを用い，界面構造の変化からポテンシャルの組み合わせを検討する．

4.2.1

シミュレーション条件

60 × 60 × 7 の単位格子からなる Si 基板上に，40 × 40 × 20 格子の GaAs を積層し，実験観察結果(14)–(18) に従い界面となす角度が 55◦ になるようにカットしたものを解析対象とする (図 4.9)．原子数は 321571 個である．Si の格子定数 laSi は 0.543 nm に固定し，堆積する GaAs の初期格子長さ laGaAs は 0.543, 0.56533, 0.58766 nm の 3 種類を考慮した．laGaAs = 0.543 nm とした場合，laGaAs = laSi であり格子ミスフィットが零で整合界面となる．ただし，0 K における GaAs の格子定数は laGaAs = 0.56533 nm である
ので，GaAs 層は圧縮状態となる．laGaAs = 0.5633 nm としたときは Si，GaAs 層とも無ひずみであるが界面は非整合となる．初期構造において GaAs 層が圧縮状態となる

laGaAs = 0.543 nm の場合 (GaAs 層の初期ひずみは −0.0395) と GaAs 層が引張状態となる場合を比較するため，GaAs 層の初期ひずみが 0.0395 となる laGaAs = 0.58766 nmの場合についても検討した．セルの x, y 方向には周期境界を適用し，z 方向は自由表
面の境界条件でまず 10000 fs の初期緩和計算を，前節で示したパラメータの取り扱い(1)〜(4) でそれぞれ行った．このとき，x, y 方向の Si 基板の応力が零となるようにセル寸法を制御している．また，GaAs に比べて Si 基板は変形しにくいと考えられるので，本研究では Si 基板の板厚方向は半無限体として，Si 層の下端面原子は端面法線方向の運動を拘束する境界条件を適用した．また，熱揺動の影響を排除するため温度は10 K とし，速度スケーリングにより制御した．

√
Si のミスフィット転位の可視化については，Si のカットオフ距離 rc = 1.2×( 3 laSi /4) =4.2 界面原子構造の評価

55

0.282 nm を基準として，カットオフ距離内に他の Si 原子が幾つ存在しているか (配位数) を算出し，その数によって判別した．通常，結晶内部では 4 つ，表面では 2 つの配

s)
GaA
(
s
c
tti
40 la
55°

60 la

ttics

40 l

attic

s (G

tics
0 lat

(Si)

6

z [001]y [010]x [100]

Fig.4.9 Hut cluster model

aAs
)

(Si)

7 lattics (Si)

位数を持つので，それ以外の配位数を持つ Si 原子を欠陥原子として可視化した．

4.2 界面原子構造の評価

4.2.2

56

シミュレーション結果および考察

緩和計算後の Si 基板の原子構造を，欠陥原子を濃く着色して図 4.10 に示す．Si/GaAs
界面において GaAs 側の欠陥原子は Si 側のそれと同じ位置に存在する．また，GaAs層の堆積により GaAs 底面のエッジ部分にひずみが生じるため，それを観察するには

GaAs 底面積よりも広範囲を観察する必要がある．そのため図 4.10 には GaAs 層を表示せずに Si 基板側を示している．また，この図では異なる初期格子長さ (laGaAs ) を用いた結果を並べて示している．Si の基板に T2 を用いたパラメータ (1) および (2) の図(a) および (b) では，基板上に欠陥原子の並びが見えるが，T3 を用いた図 (c)，(d) では欠陥はほとんど認められない．初期格子長さによる違い (初期の基板と GaAs クラス
ター間のひずみ) はポテンシャル (1) と (2) で似た傾向を示している．整合させた状態(laGaAs = 0.543 nm) から緩和するとクラスターの縁に欠陥原子が集まり，Si/GaAs 界面には欠陥原子は存在しない．一方，laGaAs = 0.56533 nm や laGaAs = 0.58766 nm(GaAsに引張ひずみ) で緩和を行った場合，界面上に平行に並ぶ転位線が発生する．図 (b) に示すようにポテンシャル (2) を用いた場合に転位線の一部が二重となっているが，図(a) のポテンシャル (1) では一本となっている．
この転位線を詳細に検討するため，図 4.10(a) の中央の図において，白線の四角で囲った領域を拡大し，2 方向から見て図 4.11 に示す．図中に白い矢印で示したのが欠陥と判定された原子である．図 4.11(a) では，ダングリングボンドを有する原子群 (Si 基
板の上端面原子) が黒く着色した原子側に引き寄せられるように配列している．また，

90◦ 回転させた図 (b) では，欠陥と判定された原子は基板内側にわずかに沈んでいる．
これは図 (a) でも破線で比較したように確認できる．次に GaAs 層を除かずに図 4.11と同じ方向，同じスケールで図 4.12 に示す．白い矢印で示しているのが図 4.11 において示した Si の欠陥原子である．図 4.12(a) の中央部分に示すように，GaAs 側に七員環が，Si 側に五員環が生じている．このため，GaAs の最下端にある Ga 原子は，Si 基板の上端面原子と同じく，図 (a) の左右方向にわずかにひずんでいるが，図 (b) を見ると

Si 基板側への沈み込みはなく，また [110] 方向 (左右方向) にもズレはない (図 4.12(b))．
次に，図 4.10(b) の中央の図において，白線の四角で囲った領域を拡大し，2 方向から

4.2 界面原子構造の評価

[001]
[010]

[100]

laGaAs= 0.543 nm

laGaAs= 0.56533 nm

laGaAs= 0.58766 nm(a) potential (1)laGaAs= 0.543 nm

la GaAs= 0.56533 nm

laGaAs= 0.58766 nm(b) potential (2)laGaAs= 0.543 nm

laGaAs= 0.56533 nm

laGaAs= 0.58766 nm(c) potential (3)laGaAs= 0.543 nm

la GaAs= 0.56533 nm

laGaAs= 0.58766 nm(d) potential (4)Fig.4.10 Defect structure on Si substrate under GaAs hut cluster

57

4.2 界面原子構造の評価

[001]

58

[001]
[110]

[110]

[110]
[110]

(a) Rotated view 1(b) Rotated view 2

Fig.4.11 Magniﬁed view of defect Si atoms and substrate indicate by white
rectangular in Fig.4.10(a)

[001]

: Si

[001]

: Ga
[110]
[110]

[110]
[110]

(a) Rotated view 1: As(b) Rotated view 2

Fig.4.12 Magniﬁed view of atoms indicate by white rectangular in Fig.4.10(a)見て図 4.13 に示す．この図には GaAs 層を除かずに示している．白い矢印および黒い矢印で示したのが欠陥と判定された原子である．白い矢印で示した原子はポテンシャ
ル (1) を用いた場合の図 4.12(a) と同じ配列を見せているが，図 4.12(a) と異なり，黒い矢印で示した原子は図の中央に引き寄せられていない．図 4.10(b) の中央の図において，白線の四角で囲った領域から [11̄0] 方向にわずかにずらした領域を拡大し，図 4.14に示す．図 4.14 において白色および黒色の矢印で示した原子は図 4.13 において示した原子と同じものである．黒い矢印で示した原子がその右側の原子に引き寄せられている．そのため，図 4.13 および図 4.14 に丸印で示した部分に原子が粗となる領域ができる．なお，図 4.10 で観察したこれらの転位線の方向が同じなのは，Si 基板上面のダングリングボンドによる．図 4.15 に示したように，今回 Si 基板上面は，図 (a) の濃く着

4.2 界面原子構造の評価

[001]

59

: Si

[001]

: Ga
[110]

[110]

[110]

[110]

(a) Rotated view 1: As(b) Rotated view 2

Fig.4.13 Magniﬁed view of atoms indicate by white rectangular in Fig.4.10(b)

[001]
[110]
[110]

Fig.4.14 Magniﬁed view of atoms indicate by black rectangular in Fig.4.10(b)色した原子に対応するので，ダングリングボンドは [11̄1] および [1̄11] 方向である．図(b) の黒色原子位置を Si 基板表面とすれば，ダングリングボンドは [111] および [1̄1̄1]方向であり，観察される欠陥も図 4.10 の転位線と直交する方向になる．
次に，T3 を用いた場合に Si/GaAs 界面上に欠陥原子が生じなかったことについて
詳細に検討する．図 4.10(c) および (d) の中央の図において，白線の四角で囲った領域を拡大し，2 方向から見て図 4.16 および図 4.17 に示す．この図には GaAs 層を除かずに示している．図 4.16(a) および図 4.17(a) に示すように，Si–Ga 結合が図の中央では[001] 方向に平行となり (白い矢印で示した結合)，図の左右ではわずかに傾いている (黒い矢印で示した結合)．図 4.10(c) の中央の図において，白線の四角で囲った領域から[11̄0] 方向にずらした黒線の四角で囲った領域を拡大し，図 4.18 に示す．図中に破線で示すように，GaAs の最下端にある Ga 原子および Si 基板の上端面原子は図の右左で

4.2 界面原子構造の評価

la

la

[001]
[010]

60

[001]
[100]

[010]

la

[100]

la

(a)

la

la

(b)

Fig.4.15 Unit lattice of Si diamond structure

[001]

: Si

[001]

: Ga
[110]
[110]

[110]
[110]

(a) Rotated view 1: As(b) Rotated view 2

Fig.4.16 Magniﬁed view of atoms indicate by white rectangular in Fig.4.10(c)

[001]

: Si

[001]

: Ga
[110]
[110]

[110]
[110]

(a) Rotated view 1: As(b) Rotated view 2

Fig.4.17 Magniﬁed view of atoms indicate by white rectangular in Fig.4.10(d)4.2 界面原子構造の評価

61

[001]
[110]
[110]

Fig.4.18 Magniﬁed view of atoms indicate by black rectangular in Fig.4.10(c)高さがわずかに異なる．しかし，Si 部，GaAs 部のみを観察すると，原子間距離がわずかに変化しているが，結合角は初期のそれとほとんど変化していない (図 4.16 および図 4.18)．そのため，各原子の配位数による欠陥判別では転位を見出せない．しかしながら，図 4.18 から分かるように界面上下の Si/GaAs 格子を通過して見ると界面に格子不整合としての転位は存在する．

T2 は表面構造をよく再現するようにパラメータフィッティングされているため，エネルギーは正しく評価されないが，ボンドの再構成が再現できる．一方，T3 は配位数による結合角度への変化が弱いため，バルクの結合角を保持し，界面が乱れない．

4.3 引張シミュレーション

4.3

62

引張シミュレーション

前節ではハットクラスターを想定したシミュレーションを行ったが，本節では前章の Ni/Ni3 Al 界面の解析と同様に，無限平面とした Si/GaAs 界面を対象として，界面法線方向に引張りを与えるシミュレーションを各ポテンシャルセットで行う．

4.3.1

シミュレーション条件

Si 層の原子面数に対して一原子面少ない GaAs 層を接合して半整合界面を作成した(図 4.19(a))．Si 層は 25 × 25 × 12 格子，GaAs は 24 × 24 × 12 格子で作成し，原子数は 115296 個である．このセルの x, y 方向に周期境界を適応し，z 方向は自由表面の境界条件で 10000 fs の初期緩和計算を行った．このとき，x, y 方向の応力が零となるよ
うにセル寸法を制御している．また，熱揺動の影響を排除するため温度は 10 K とし，
速度スケーリングにより制御した．ポテンシャルパラメータに関しては前節と同じく，

4.1.2 項に示したポテンシャル (1)〜(4) の 4 種類を考慮する．ポテンシャル (1) および(3) を用いたモデルにおいて初期緩和計算後に欠陥原子と判定した原子を濃く着色して Si 層とともに図 4.19(b) および (c) に示す．前節と同じく，ポテンシャル (1) を用いた場合は Si/GaAs 界面上に平行に並ぶ転位線が生じ (図 (b))，ポテンシャル (3) を用いた場合は界面上に欠陥原子が生じていない (図 (c))．緩和計算後は，x, y 方向の応力が零になるように制御しながら，z 方向に微小ひずみ増分 dϵ = 1.0 × 10−6 [ 1 / fs ] を毎ス
テップ与えて引張るシミュレーションを行った．

4.3 引張シミュレーション

24 l
att
(13. ices (G
6 nm
aAs
)
)

)

aAs
s (G
e
c
i
t
t
4 la

12 lattices(6.78 nm)2

12 lattices(6.52 nm)

GaAs

Si

z [001]y [010]25 l
att(13. ices (Si
6 nm)
x [100]
)

es
attic
25 l

(Si)

(a) Simulation model(b) Defect structure on Si/GaAs
interface by potential (1)(c) Defect structure on Si/GaAs
interface by potential (3)Fig.4.19 Simulation cell and defect structures after relaxation

634.3 引張シミュレーション

4.3.2

64

シミュレーション結果および考察

4 種類のパラメータセットにおいて，それぞれ引張りを行った際の応力 σz とひずみ
ϵ の関係を図 4.20 に示す．Si のパラメータに T2 を用いた場合 (図 (a)，(b)) よりも T3を用いた場合 (図 (c)，(d)) の方が最大応力を生じるひずみが小さい．ピーク応力後，いずれも応力急減して破断した．破断した原子構造を図 4.21 に示す．図 (a) および (b)は T2 を用いたときの ϵ = 0.135 における原子状構造，図 (c) および (d) は T3 を用いたときの ϵ = 0.115 および ϵ = 0.110 における原子構造を示している．Si のパラメータに T2 を用いた場合には，図 (a) および図 (b) に示すように，界面から {111} 面に沿ったき裂が生じている．その様子を詳しく検討するため，ポテンシャル (1) を用いた場
合について，Si 基板および GaAs の欠陥原子のみを図 4.22 に示す．図で薄く着色した原子が GaAs の欠陥である．図には ϵ = 0.120〜0.130 における様子を並べて示している．ϵ = 0.120 では界面上に格子状の転位線が見られる．その後，ϵ = 0.125 において，転位線の交点から GaAs 層側に欠陥原子が多数発生し (図の丸印)，ϵ = 0.130 になると {111} 面上に成長している．この欠陥はき裂部分でボンドが切れた原子である．周期境界条件下では，すべりが交差する部分からへき開が発生することが報告されている(42) ．本シミュレーションでも，界面転位の交点からへき開を生じた後，(111) 面に成長したものと理解できる．一方，Si のパラメータに T3 を用いた場合には，界面を起点とせずに GaAs 層内の様々な所で {111} 面でのへき開を生じている．

4.3 引張シミュレーション

64

ε = 0.135
2

0

0.05

0.1

Stress σ zz [GPa]Stress σ zz [GPa]

6

0.15

ε = 0.135

4

2

0

0.05

(a) potential (1)(b) potential (2)6Stress σ zz [GPa]6Stress σ zz [GPa]0.15Strain εStrain ε4ε = 0.115

2

0

0.1

0.05

0.1

Strain ε(c) potential (3)

0.15

4

ε = 0.110
2

0

0.05

0.1

Strain ε(d) potential (4)Fig.4.20 Relationship between stress and strain under tension

0.15

65

4.3 引張シミュレーション(a) potential (1), ϵ = 0.135(b) potential (2), ϵ = 0.135(c) potential (3), ϵ = 0.115(d) potential (4), ϵ = 0.110

Fig.4.21 Atomic structure in yielding by tension simulation

ε = 0.120

ε = 0.125

ε = 0.130

Fig.4.22 Si substrate and defect atoms of GaAs (potential (1))

66

4.4 結言

4.4

67

結言

Si には Tersoﬀ 型ポテンシャル，GaAs には Jeﬀerson が Tersoﬀ 型ポテンシャルにフィッ
ティングしたパラメータがあるものの，Si/GaAs 界面の原子シミュレーションを行う
とき，界面におけるポテンシャルパラメータの取り扱いが不明である．本章では，まず膨張圧縮に対する Si，GaAs 単相のエネルギー変化を第一原理計算と比較した後，Si，

GaAs 単位格子を積層したスーパーセルで界面エネルギーを求めた．Si/GaAs 界面での Tersoﬀ 型ポテンシャルの問題は，原子 i を中心として，他の原子 j, k との結合角を考慮するポテンシャルパラメータが，Tersoﬀ の Si では原子 i の種類だけに依存するのに対し，Jeﬀerson の GaAs では原子 i と j の組に陽に依存する点である．そこで，[ 1 ]結合角をなす原子 i, j, k が Ga と As のみのときは Jeﬀerson の Ga–As にフィッティン
グされたパラメータを用い，1 つでも Si 原子が入った場合は原子 i の単相 (Si，Ga，Asのいずれか) のパラメータとする，[ 2 ] 基本のボンドを構成する原子 i, j が Ga と Asであるときは，原子 k が Si でも Jeﬀerson の Ga–As パラメータを用い，それ以外は原子 i に依存させる，の 2 通りを考慮して界面エネルギーを算出し，第一原理の計算結果と比較した．なお，Si のポテンシャルパラメータには T2 および T3 と呼ばれる 2 種類のパラメータが提案されているので，計 4 通りのポテンシャルセットで比較した．その後，Si 基板上のハットクラスターや無限平面とした Si/GaAs 界面の MD シミュレー
ションを行い，4 通りのポテンシャルパラメータで観察される Si/GaAs 界面の構造や強度などを検討した．得られた結果を要約して以下に示す．

I. 界面エネルギー(1) Si と GaAs が混在する部分の結合角について，ポテンシャル [ 1 ] の取り扱いをすると，T2，T3 いずれも圧縮ひずみ下の界面エネルギーが非常に大きくなり第一原理の計算結果から大きく逸脱する．T2 で結合角を [ 2 ] として取り扱ったときのひずみに対する界面エネルギーの変化は第一原理の傾向に最も近いが，界面エネルギーは負である．したがって，界面エネルギーの点からは，Si に T3 を用
い，結合角には [ 2 ] の取り扱いをするのが良いという結果になった．

4.4 結言

68

II. GaAs ハットクラスター下の界面原子構造(2) 結合角の取り扱い [ 1 ]，[ 2 ] に関わらず，Si に T3 を用いた場合は Si 基板の結晶構造には乱れを生じない．一方，T2 では [ 1 ]，[ 2 ] いずれも Si 基板に欠陥を生じた．GaAs の格子長さを Si と同じに圧縮して設置したハットクラスターでは，
縁に欠陥原子が集まりクラスター下部の Si/GaAs 界面上には欠陥原子は発生しない．無ひずみおよび GaAs を膨張させて Si 基板に乗せたクラスターでは，Si 基
板に五員環，GaAs 底面に七員環が生じ，界面に一定間隔で並ぶ転位線が認められた．

III. 界面の剥離シミュレーション(3) Si に T2 を用いた場合は，結合角の取り扱いに関わらず，界面から {111} 面に沿ったき裂が生じて破断した．これはすべり交差部からへき開を生じ，それが {111}
面方向に成長したためである．一方，T3 を用いた場合は，界面から離れた GaAs
層内に {111} 面のへき開を生じた．

第5章
結論

本研究は，異種材料界面に形成される転位をテーマに，分子動力学シミュレーションを行った成果をまとめたものである．その対象には，Ni 基単結晶超合金 γ/γ ′ 界面，ならびに，半導体デバイスの Si/GaAs 界面を取り上げた．特に，Si/GaAs 界面のシミュ
レーションでは，Tersoﬀ 型ポテンシャルの界面での取り扱いを明確にすべく，第一原
理計算との比較も行っている．以下に，得られた結果を総括する．
第 2 章では，解析手法の概要を述べた．はじめに，分子動力学方の基礎方程式を示
し，原子間相互作用の評価に用いた原子埋め込み法ポテンシャルおよび Tersoﬀ 型ポテ
ンシャルの概要を説明した．次に，分子動力学計算における温度制御方法と大規模分
子動力学計算に用いられる領域分割による高速化について述べた．
第 3 章では，Ni 基単結晶超合金 γ/γ ′ 界面上のミスフィット転位と，単相中の刃状およびらせん転位の転位芯レベルの切り合い，ならびに，外力下でのミスフィット転位の運動について検討した．ミスフィット転位と単相中を進行する転位の相互作用は，基本
的にはバーガースベクトルの可算則 (転位論) に従った反応をする．しかしながら，ミスフィット転位と，同符号側に余剰原子面を持つ刃状転位が衝突すると，反力により⟨110⟩ 方向のミスフィット転位を交線とする 2 つのすべり面上に，SFT(stacking fault
tetrahedron) のような三角屋根形状の積層欠陥の生成が見られた．この 2 つのすべり面
は，進行してきた刃状転位のそれとは等しくないが，片方の面上では完全転位のバー
ガースベクトルが刃状転位のそれと一致して，4 つの部分転位の交点がすべり面の交線上をスライドすることで運動するメカニズムを見出した．また，γ/γ ′ 界面に [100] 方

69

70
向または [110] 方向にせん断を与えるシミュレーションを行い，[100] 方向へのせん断
では，界面転位網は界面上に多数の点欠陥を残しながら全体的にすべること，[110] 方
向へのせん断では，転位網の横糸部分のみ界面をすべること，などを明らかにした．

Si，GaAs それぞれには Tersoﬀ 型ポテンシャルが提案されているが，先述のように
Si/GaAs 界面でのポテンシャルパラメータの取り扱いは不明である．そこで第 4 章では，まず Si と GaAs の単位格子を接合したスーパーセルについて，VASP を用いた第一原理計算により界面エネルギーを検証した．その後，Si 基板上に GaAs ハットクラスターを種々の格子長さで設置し，Si の 2 種類の Tersoﬀ 型ポテンシャルと，2 種類の

Si/GaAs 界面の取り扱い方の計 4 通りのポテンシャルパラメータの組み合わせで構造緩和を行った．界面での取り扱いによる差は大きくなく，Si 部に表面構造を再現する

T2 を用いるか，弾性係数を良好に再現する T3 を用いるかが界面の欠陥構造形成には重要となった．T2 ではクラスターの縁または下の界面に，Si の再構成による欠陥を生じる．一方，T3 を用いた場合は Si，GaAs 層それぞれは単位結晶でのダイヤモンド構
造から大きく変化することはない．最後に，無限平面とした Si/GaAs 界面に，界面法
線方向へ引張を与えるシミュレーションを行った．Si に T2 を用いた場合は，Si 基板
上の初期欠陥を起点として，界面から {111} 面に沿った亀裂が生じて破断した．一方，

T3 を用いた場合は，界面から離れた GaAs 層内に {111} 面のへき開を生じた．
以上，金属の半整合界面と，半導体の異種材料界面について分子動力学シミュレー
ションによる様々な検討を行った．Ni 基超合金について得た知見は，離散転位動力学
法 (Discrete Dislocation Dynamics; DDD)(6), (7) によるスケールアップなどで，Ni 基超合金の力学的階層性の解明へ繋がることが期待される．一方，半導体 Si/GaAs 界面のシミュレーションでは，原子間ポテンシャルを仮定したシミュレーションの限界も露呈した．Si/GaAs 界面の正確なシミュレーションを行うためには，Si–Ga，Si–As の結合パラメータについて，第一原理計算によるさらなる検証が必要と考える．

参考文献(1) Pollock, T. M. and Argon, A. S., Creep Resistance of CMSX–3 Nickel BaseSuperalloy Single Crystals, Acta Metall. Mater., 40–1 (1992), 1–30.
(2) Zhang, J. X., Murakumo, T., Harada, H. and Koizumi, Y., Dependence of creep
strength on the interfacial dislocation in a fourth generation SC superalloy TMS–
138, Scripta Materialia, 48 (2003), 287–293.
(3) Zhang, J. X., Murakumo, T., Koizumi, Y. and Harada, H., The inﬂuence ofinterfacial dislocation arrangements in a fourth generation single crystal TMS–
138 superalloy on creep properties, Journal of Materials Science, 38 (2003),
4883–4888.
(4) Zhang, J. X., Murakumo, T., Koizumi, Y., Kobayashi, T., Harada, H. andMasaki, S., Interfacial dislocation networks strengthening a fourth–generation
single-crystal TMS-138 superalloy, Metallurgical and Materials Transactions A,
33 (2003), 3741–3746.
(5) Zang, J. X., Koizumi, Y., Kobayashi, T., Murakami, T. and Harada, H., Mechanical Behavior–Communication : Strengthening by γ/γ ′ Interfacial Dislocation Networks in TMS–162 – Toward a Fifth–Generation Single–Crystal Superalloy, Metallurgical and Materials Transactions – Series A, 35–6 (2004), 1911–
1914.
(6) Zibib, H. M., Rhee, M. and Hirth, J. P., On Plastic Deformation and The
Dynamics of 3D Dislocations, Int. J. Mech. Sci., 40–2–3 (1998), 113–127.
(7) Van der Giessen, E. and Needleman, A., Discrete Dislocation Plasticity, Modeling Simulation Mater. Sci. Eng., 3 (1995), 689–735.

71

参考文献

72

(8) 相原智泰, 正朋祥, 川添良幸, 傾斜構造を持つ Ni/Ni3 Al 界面の弾性挙動の分子動力学シミュレーション, 日本金属学会誌, 62–11 (1998), 978–985.

(9) 屋代如月, 内藤正登, 冨田佳宏, Ni 基超合金格子状析出相の頂点近傍における転位の発生および運動の分子動力学解析, 日本機械学会論文集, A 編, 68–668 (2002),

573–579.
(10) 屋代如月, 内藤正登, 冨田佳宏, Ni 基超合金格子状析出相に堆積する刃状転位の分子動力学シミュレーション, 日本機械学会論文集, A 編, 69–677 (2003), 210–217.

(11) Zhu, T. and Wang, C., Misﬁt dislocation networks in the γ/γ ′ phase interface
of a Ni–based single–crystal superalloy : Molecular dynamics simulations, Phys.
Rev. B, 72 (2005), 014111.
(12) 屋代如月, Pangestu, J. R., 冨田佳宏, Ni/Ni3 Al 界面ミスフィット転位の構造および運動 : 分子動力学による検討, 材料, 56–5 (2007), 439–444.

(13) Yashiro, K., Pangestu, J. R. and Tomita, Y., Molecular Dynamics Study of Interfacial Dislocation Network at Gamma/Gamma–Prime Interface in Ni–Based
Superalloys，Third International Conference on Multiscale Materials Modeling
MMM2006 Freiburg, (2006), 582–584.
(14) Usui, H., Yasuda, H. and Mori, H., Morphology and lattice coherency in GaAs
nanocrystals grown on Si(100) surface, Appl. Phys. Lett. 89 (2006), 173127.
(15) Usui, H., Ishiji, K., Yasuda, H. and Mori, H., Growth process and nanostructureof crystalline GaAs on Si(110) surface prepared by molecular beam epitaxy,
Journal of Crystal Growth, 289–1 (2006), 76–80.
(16) Usui, H., Yasuda, H. and Mori, H., Crystal shape of GaAs nanocrystals deposited on Si(100) by molecular beam epitaxy, Thin Solid Films, 516–9 (2008),
2487–2490.

参考文献

73

(17) Usui, H., Yasuda, H. and Mori, H., Eﬀect of substrate temperature on growth
process of GaAs on Si(110) vicinal surface studied by reﬂection high–energy
electron diﬀraction, Physica E : Low–dimensional Systems and Nanostructures,
40–4 (2008), 836–842.
(18) 福田幸正, Si(100) 表面上にヘテロ成長した InAs ナノ結晶の形状異方性, 神戸大学大学院 工学研究科 機械工学専攻 修士論文, (2008).

(19) 森田達夫, 半導体薄膜形成のためのリモートラインソースプラズマ CVD 装置の研究開発及びそれに基づく起業実践, 高知工科大学大学院 工学研究科 基盤工学
専攻 起業家コース 博士 (工学) 学位論文, (2001).

(20) 守谷浩志, 岩崎富生, 保川彰夫, 三浦英生, 分子動力学を用いた Si 薄膜真性応力と微細構造の解析, 日本機械学会論文集, A 編, 63–613 (1997), 1999–2005.

(21) 小山紀久, 岡島康, 太田英二, 武田京三郎, 白石賢二, 山口浩司, 伊藤智徳, 大野隆
央, 自己組織化量子ドット形成におけるミスフィット転位の影響, 第 60 回応用物
理学会学術講演会, (1999).

(22) 石井晃, 三宅大助, 逢坂豪, 第一原理計算による GaN(0001) のエピタキシャル成長の研究, 日本結晶成長学会, 29–2 (2002), 82.

(23) 伊藤 智徳, 白石 賢二, 田口 明仁, GaAs 表面上吸着過程に関する量子論的シミュ
レーション, 日本結晶成長学会誌, 26–2 (1999), 49–50.

(24) Daw, M. S. and Baskes, M. I., Semiempirical Quantum Mechanical Calculationof Hydrogen Embrittlement in Metals, Phys. Rev. Lett., 50–17(1983), 1285–
1288.
(25) Daw, M. S. and Baskes, M. I., Embedded–atom method : Derivation and Application to Impurities, Surfaces and Other Defects in Metals, Phys. Rev. B, 29–12(1984), 6443–6453.

参考文献

74

(26) Tersoﬀ, J., New empirical approach for the structure and energy of covalentsystems, Phys. Rev. B, 37–12 (1987), 6991.
(27) Tersoﬀ, J., Empirical interatomic potential for silicon with improved elastic
properties, Phys. Rev. B, 38–14 (1988), 9902.
(28) Tersoﬀ, J., Modeling solid–state shemistry : Interatomic potentials for multicomponent systems, Phys. Rev. B, 39–8 (1988), 5566.
(29) Voter, A. F. and Chen, S. P., Accurate Interatomic Potentials for Ni, Al andNi3 Al, Mat. Res. Soc. Symp. Proc., 82 (1987), 175–180.
(30) Chen, S. P., Voter, A. F., Albers, R. C., Boring, A. M. and Hay, P. J., Investigation of the Eﬀect of Boron on Ni3 Al Grain Boundaries by Atomistic Simulations,J. Mater. Res., 5–5 (1990), 955–970.
(31) Voter, A. F., Chapter 4 The Embedded Atom Method, John Wiley & Sons,(1994), 77–90.
(32) Foiles, S. M., Calculation of the Surface Segregation of Ni–Cu Alloys with theUse of the Embedded–Atom Method, Phys. Rev. B, 32–12 (1985), 7685–7693.
(33) Rose, J. H. Smith, J. R., Guinea, F. and Ferrante, J., Universal features of theequation of state of metals, Phys. Rev. B, 29 (1984), 2963–2969.
(34) Balamane, H., Halicioglu, T. and Tiller, W. A., Phys. Rev. B, 46 (1992), 2250.
(35) 丸山茂夫, 井上知洋, 熱流体系および固体系のミクロシミュレーションに関する
合同シンポジウム・第 5 回分子動力学シンポジウム, (2000), 11–12.

(36) Smith, R., Nucl. Instr. and Meth., Phys. Res. B, 67, 335 (1992).
(37) Sayed, M., Jeﬀerson, J. M., Walker, A. B. and Cullis, A. G., Molecular dynamics simulations of implantation damage and recovery in semiconductors, Nucl.
Instrum. Meth. Phys. Res. B, 102 (1995), 218.

参考文献

75

(38) 上田顯, コンピューターシミュレーション, (1990), 朝倉書店.
(39) Honeycutt, J. D. and Andersen, H. C., Molecular Dynamics Study of Melting
and Freezing of Small Lennard–Jones Clusters, J. Phys. Chem., 91 (1987), 4950–
4963.
(40) Vanderbilt, D., Phys. Rev. B, 4, (1990), 7892.
(41) Kresse, G. and Furthmüller, J., Phys. Rev. B, 54, (1996)，11169.
(42) 屋代如月, 北村隆行, 大谷隆一, 引張軸に対して直角方向に拘束を受けるニッケル
結晶のへき開破壊に関する原子レベル解析, 日本機械学会論文集, A 編, 63–615(1997), 2380–2386.

関連発表論文・講演論文
関連発表論文¤ Kisaragi Yashiro, Yuuhu Suzuki, Joy Rizki Pangestu and Yoshihiro Tomita,
Molecular Dynamics Study on Characteristics of Misﬁt Dislocations in Ni–Based
Superalloys, Key Engineering Materials Vols., 345–346 (2007), 951–954.
¤ 鈴木雄風, 屋代如月, 冨田佳宏, 刃状およびらせん転位とミスフィット転位の転位
芯相互作用：分子動力学による解析, 日本機械学会論文集, A 編, 73–735 (2007),

1217–1224.
講演論文¤ 鈴木雄風, 屋代如月, 冨田佳宏, 界面転位網と刃状およびらせん転位の転位芯相互作用, 第 12 回分子動力学シンポジウム, (2007.5).

¤ Kisaragi Yashiro, Yuuhu Suzuki, Joy Rizki Pangestu and Yoshihiro Tomita,
Molecular Dynamics Study on Characteristics of Misﬁt Dislocations in Ni–Based
Superalloys, 10th International Conference on the Mechanical Behavior of Materials, (2007.5).
¤ Kisaragi Yashiro, Yuuhu Suzuki, Joy Rizki Pangestu and Yoshihiro Tomita,
Molecular Dynamics Study on Core Interaction between Edge/Screw Dislocation
and Interfacial Dislocation Network, The International Conference on Structural
Analysis of Advanced Materials, (2007.9).
¤ 鈴木雄風, 屋代如月, 冨田佳宏, 分子動力学法による Si/GaAs 界面エネルギーの評価, M&M2008 材料力学カンファレンス, (2008.9).

76

関連発表論文・講演論文
機論-1/877

関連発表論文・講演論文
機論-2/878

関連発表論文・講演論文
機論-3/879

関連発表論文・講演論文
機論-4/880

関連発表論文・講演論文
機論-5/881

関連発表論文・講演論文
機論-6/882

関連発表論文・講演論文
機論-7/883

関連発表論文・講演論文
機論-8/884

関連発表論文・講演論文

MD シンポジウム-1/6

85

関連発表論文・講演論文

MD シンポジウム-2/6

86関連発表論文・講演論文

MD シンポジウム-3/6

87関連発表論文・講演論文

MD シンポジウム-4/6

88関連発表論文・講演論文

MD シンポジウム-5/6

89関連発表論文・講演論文

MD シンポジウム-6/6

90関連発表論文・講演論文

ICM10-1/4

91

関連発表論文・講演論文

ICM10-2/4

92

関連発表論文・講演論文

ICM10-3/4

93

関連発表論文・講演論文

ICM10-4/4

94

関連発表論文・講演論文

Edge/Screw-1/1

95

関連発表論文・講演論文

M&M-1/2

96

関連発表論文・講演論文

M&M-2/2

97

謝辞
本研究を遂行するにあたり，冨田佳宏教授の広い御見識と暖かい御人柄に導かれた
ことは，著者にとって無類の幸福でありました．また，教授の懇切丁寧な御指導に最大限の敬意と感謝の念を表すとともに心より御礼申し上げます．浅学非才の著者が本論文を完成させるまで研究全般に渡り直接御指導を頂いたのみならず，日頃の研究生
活においても多くの御助言・御協力頂きました屋代如月准教授に深く感謝いたします．
研究活動に対して幅広い視点からの御助言により，著者の狭い見解を広げて下さった
長谷部忠司准教授に感謝いたします．また，日々の研究活動が円滑に行えるよう環境づくりに御尽力下さった古宇田由夫技術職員に心から御礼申し上げます．
研究室に配属後，共に研究活動に取り組んだ小西正彰君，近堂将規君，竹内耕助君，
久國陽介君，廣内智之君，山田将臣君，山本雄太君，横川望君に感謝いたします．さらに同じ研究グループの中島裕介さん (現，川崎重工業株式会社)，西村正臣さん，Joy

Rizki Pangestu さんをはじめとする固体力学研究室内諸先輩方に深く感謝いたします．
最後に，大学生活において積極的に支援して頂いた家族に心から感謝の意を表します．ありがとうございました．

2009 年 2 月
鈴木98

雄風