修士論文

シリコンの塑性変形メカニズムの基礎的・定量的検討：原子弾性剛性係数による局所格子不安定性解析

指導教員：屋代 如月

藤原 正大2013 年 2 月
神戸大学大学院 工学研究科 博士課程前期課程 機械工学専攻

Master Thesis
Quantitative Approarch on Primitive Mechanismof Plastic Deformation in Silicon:
Local Lattice Stability Analysis Based on Atomic
Elastic Stiﬀness Coeﬃcients

February 2013
Division of Mechanical Engineering,
Graduate School of Engineering,
Kobe University, Kobe, Japan
MASAHIRO FUJIHARA

要約
本論文では，分子動力学法によりシリコンの単結晶バルク，薄板，ナノワイヤそして粒界積層構造，バンブー構造などについて変形シミュレーションを行い，不安定挙動(応力急減) と原子弾性剛性係数 (Bijα ) の関係について検討した．大きさを統一した小規
模な周期セルで，温度を変えた単結晶バルク，ナノワイヤ，粒界積層構造などの [001]
方向引張シミュレーションを行った．完全結晶に近い状態では系の不安定（応力急減）
前に応力上昇の頭打ち→なだらかな低下が存在し，その時に detBijα の揺らぎが増大していること，熱揺らぎや構造不均一を有する系でも detBijα のばらつきは引張により減少するが，応力ピークより前のひずみで再び増大すること，などを明らかにした．次に表面による構造不均一をより詳細に議論するため，より大きなセルで (100) および(110) 表面の薄板，ナノワイヤについて [001] 方向引張シミュレーションを行った．やはり引張により標準偏差の幅が減少するが，先に述べたようなエラーバーの拡大を生じないケースもあり，(100) 表面ナノワイヤの結果からは，応力急減直前のエラーバーのパルス的な増大はエッジ部に大きな負の detBijα が出現したためであることが示唆された．最後に異なる変形モードであるせん断変形シミュレーションを行い，ピーク直後のひずみにおいて，内部にせん断帯状の detBijα < 0 の原子群が発生し，そこにすべりを生じて [1̄1̄0] 方向に連続した detBijα < 0 の原子群を形成するのが確認された．また変形速度が遅いシミュレーションの詳細な観察結果により，応力上昇が頭打ちしゆらぎを生じる点=detBijα < 0 の原子が現れる点であること，構造変化を生じ応力がカ
タストロフ的に急減する点=負の detBijα がゆらぎながら増加し，最終的に大きな負値を生じる「不安定」点であることが示唆された．さらに Bijα の固有値を調べ，変形直前に Bijα の負の最小の固有値の固有ベクトルの対称性が失われていることを示した．

Summary
In order to elucidate the relationship between unstable stress drop and atomic elastic
stiﬀness( Bijα ), we performed various deformation simulations on bulk, thin plate, nanowire, and laminate/bamboo structure of grain boundary for silicon using molecular dynamics method. First, we have performed various tensile simulations on bulk/nanowireof Si single crystal, laminate-bulk/bamboo-nanowire with Σ5 twist grain boundary. It
revealed that there is a slight and smooth stress peak before the unstable stress drop,
and the standard deviation of det Bijα began to increase at the point, in the case ofbulk with very low temperature. Even in the systems with thermal fluctuation andstructural inhomogeneity, the standard deviation of det Bijα decreases at the initialstage of tension, but it increases again precursory before the unstable stress drop. In
order to discuss the eﬀect of structural inhomogeneity by surface more detail, we thenperformed [001] tensile simulations on thin plate and nanowire for (100) and (110)
surfaces with larger cell. As same as results mentioned above, the standard deviationof det Bijα decrease by tension, while the precusory fluctuation increase can not befound in some case. From the results of the (100) surface nanowire, it is suggested
that pulse-like increase in the standard deviation of det Bijα just before the unstablestress drop is attributed to the emergence of large negative value of det Bijα on theedge of nanowire. Finally, we performed shear simulation on a Si single crystal as
the diﬀerent deformation mode, we found that shear-band like domains of det Bijα < 0atoms emerge just after the stress-strain peak, then the bands change its orientation to[1̄1̄0] direction due to the slip at the band. Detailed observation of the simulation withslower deformation speed revealed that the stress begin to fluctuate at the point whendet Bijα < 0 atoms emerge, and the system shows catastrophic stress drop at the pointof ”instability”, where det Bijα < 0 emerges in the fluctuation of emerge/disappear ofdet Bijα < 0 domain. Further, we examined the eigenvalue, and demonstrated that thesymmetry of the eigenvector smallest negative eigenvalue was lost at the point justbefore deformation.

目次
第 1 章 緒言1

第 2 章 解析手法

42.1

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

4

2.2

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

5

2.3

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

6

2.4

Tersoﬀ ポテンシャルにおける力の導出 . . . . . . . . . . . . . . . . . .

82.5

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

13

2.6

格子不安定解析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

2.6.1

Born の安定基準 . . . . . . . . . . . . . . . . . . . . . . . . . .

14

2.6.2

ブラベー格子を用いた安定解析 . . . . . . . . . . . . . . . . . .

15

2.6.3

弾性剛性係数による安定性解析 . . . . . . . . . . . . . . . . . .

15

第 3 章 小規模セルによる構造不均一性の包括的検討

3.1

3.2

3.3

18

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

3.1.1

温度の効果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

3.1.2

構造不均一性の効果 . . . . . . . . . . . . . . . . . . . . . . . .

19

解析結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

3.2.1

応力−ひずみ曲線 . . . . . . . . . . . . . . . . . . . . . . . . . .

21

3.2.2

原子弾性剛性係数のばらつき . . . . . . . . . . . . . . . . . . .

22

3.2.3

不安定開始点の特定 . . . . . . . . . . . . . . . . . . . . . . . .

23

表面・粒界の解析結果および考察 . . . . . . . . . . . . . . . . . . . . .

27

3.3.1

表面の効果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.3.2

粒界の効果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

i

目次
第 4 章 薄板・ナノワイヤによる表面不均一性の議論

ii
35

4.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

4.2

解析結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

4.2.1

(100) 表面薄板の解析 . . . . . . . . . . . . . . . . . . . . . . . .

37

4.2.2

(110) 表面薄板の解析 . . . . . . . . . . . . . . . . . . . . . . . .

41

4.2.3

(100)-(010) ナノワイヤの解析 . . . . . . . . . . . . . . . . . . .

44

4.2.4

(110)-(1̄10) ナノワイヤの解析 . . . . . . . . . . . . . . . . . . .

48

第 5 章 せん断変形下の不安定挙動の検討

52

5.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52

5.2

解析結果および考察. . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

5.2.1

温度の効果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

5.2.2

ひずみ速度の効果 . . . . . . . . . . . . . . . . . . . . . . . . . .

58

固有値解析による不安定挙動の検討 . . . . . . . . . . . . . . . . . . . .

63

5.3

第 6 章 結論

68

参考文献

70

関連発表論文・講演論文

73

謝辞

81

第1章
緒言

シリコンは IC や LSI 等のエレクトロニクス機器の基礎となる重要な材料である．その微細加工技術を背景に，シリコンは MEMS 分野への応用も多い．MEMS の代表格であるセンサは小さい方が望ましく，微小化することで，軽量化、省スペース，高速化は当然ながら，物理量や化学量を測る際，小さいほど測定対象に影響を与えにくく，
ピンポイントで測定点の正確な情報を感度良く得ることができるなど様々なメリットがある．このような動機からも，さらなる微小化を目指したマイクロ・ナノスケール
での変形特性の研究が活発に行われている．
シリコンのマイクロ・ナノスケールでの力学特性を実験的に評価した研究には，原子間力顕微鏡を利用した単結晶シリコンナノワイヤの曲げ試験や (1) ，ナノインデンテー
ション周囲の構造がアモルファスに構造変化を観察した例 (2) などがある．
近年の計算機能力の向上を背景に，原子シミュレーションによるシリコンの変形・破壊現象の解明も盛んに行われている．Jing らは分子動力学法を用いた単結晶シリコ
ンナノワイヤ (SiNWs) の単軸引張，単軸圧縮シミュレーションを行っている (3), (4) ．単
軸引張シミュレーションでは臨界荷重は温度の増加またひずみ率の増加で小さくなり，
ナノワイヤの直径の増加によって大きくなること，軸方向の結晶方位が異なる単結晶
ナノワイヤの単軸圧縮シミュレーションでは結晶方位によって変形挙動が異なること，などを報告している．Park らはアモルファスシリコンの様々な温度で結晶化過程をシ
ミュレーション (5) するとともに，シリコンナノカンチレバーの曲げと水平振動によって弾性特性のサイズ効果を検討している (6) ．泉らは微視的スケールにおいては，アモ1ルファス構造に起因した空間的なばらつきが重要になることに着目し，分子動力学法
によりアモルファスシリコンの表面エネルギー・表面応力を求めている (7) ．また香山らは Si のねじり粒界の原子・電子構造・界面エネルギー解明を試み，対称傾角粒界に比べてねじり粒界が一般に高エネルギーであることを示し，多結晶 Si 中でねじり粒界がほとんど観察されないことの理由をよく説明した (8) ．なお，対象とした原子数が少ないため直接変形・破壊を扱ったものではないが，Si の基礎的物性の評価として第一原理計算による検討も広くなされており，単結晶シリコンの理想せん断強度 (9), (10) や
理想引張強度等 (11) が多数報告されている．
本研究ではバルク，表面薄板，ナノワイヤ，粒界積層構造など様々な形状のシリコンについて単軸引張やせん断変形シミュレーションを分子動力学法を用いて行い，不安定挙動を観察するとともに，局所格子不安定解析を行い，不安定挙動と原子弾性剛性係数の関係について検討する．局所格子不安定解析とは各原子位置におけるエネル
ギーの空間勾配に相当する物理量である原子弾性剛性係数 (Bijα ) の正値性に着目し，局所の安定性を判別する手法であり，これまでの研究により detBijα の正負で単結晶表面からの転位発生 (12), (14) や微視的へき開破壊の発生 (13) などが評価できることが示されている．
第 2 章では解析手法について述べる．はじめに，分子動力学法の基礎式を示し，本研究で原子間相互作用の評価に用いた Tersoﬀ 型 (15), (16) ポテンシャルの概要について説明する．さらに各原子の安定性を評価するために用いた局所格子不安定解析について概説する．
第 3 章では Tersoﬀ 型ポテンシャル (15), (16) を用いた単結晶シリコンの [001] 引張変形
下の不安定挙動（応力急減）を分子動力学シミュレーションにより観察するとともに，そのときの各原子の原子弾性剛性係数について詳細に検討する．熱揺動や系の不均一性を考慮した包括的な議論とするため，小規模な周期セルによる単結晶バルクならびに粒界積層構造，横方向自由境界による単結晶ナノワイヤならびに粒界バンブー構造など様々な系について検討する．
第 4 章では，3 章のナノワイヤの結果は，表面と表面の重畳したエッジ部分の効果が大きいと考えられるため，セルの寸法を大きくして原子数を増やし，ナノワイヤだけ2

でなく薄板状モデルについて同様の [001] 引張シミュレーションを行い，表面を有する系の不安定挙動と原子弾性剛性係数の関係について詳細に検討する．また対象とする表面も (100) だけでなく (110) 表面についても検討する．
第 5 章では，異なる変形モードであるせん断変形の際の不安定挙動を観察し，そのときの各原子の原子弾性剛性係数について詳細に検討する．さらにこれまでの解析では各原子の Bijα の行列式である detBijα に着目して解析を行ってきたが，Bijα の固有値についても検討し，不安定挙動と Bijα の関係をより詳細に検討する．
最後に，第 6 章で本研究の総括を述べる．

3

第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 )
2
3
dt
2 dt
3! dt

(2.3)

となる．両式の和をとり式 (2.1) を代入すると

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

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

時刻 t + ∆t における任意の原子 i の位置を求めることができる．

4

(2.4)

2.2

原子間ポテンシャル

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

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

MD 法において，原子 i に作用する力 F i は系のエネルギー Etot の空間微分によって求めるため (式 (2.2))，系のポテンシャルエネルギー Etot をいかに精度よく評価するかが重要となる．(3) の第一原理分子動力学法は，計算量が極めて膨大になるため，変形・破壊のような多数の原子の動的挙動への直接的な適用は困難である．そこで，原子間相互作用を簡略評価する (1)，(2) の原子間ポテンシャルが通常用いられる．

5

2.3

Tersoﬀ 型ポテンシャル

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

E=

i

fc (rij )[fR (rij ) + bij fA (rij )]

(2.5)

j>i

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

1
E = [fR (rij ) + (bij + bji )fA (rij )]
2

(2.6)

fR (rij ) = A exp(−λ1 rij )fc (rij )(2.7)

fA (rij ) = −B exp(−λ2 rij )fc (rij )(2.8)

カットオフ関数は

fc (rij ) =



 1,

[

 1






2
0,

−

1
π(r − R)
sin
2
2D

]

r <R−DR−D <r <R+D

(2.9)

r >R+D

ここで fR (rij ) は斥力を表わす項であり， 12 (bij + bji ) は引力を表わす項である．またfc (rij ) はカットオフを表わす項である．bij は，原子 i，原子 j 以外に原子 k も含めた
3 個の原子によって定められる項であり，次式で表わされる．

bij = (1 + β n ζijn )− 2n1

6(2.10)

ζij =∑[

fc (rik )g(θ) exp λ33 (rij − rik )3](2.11)

k̸=i,jg(θ) = 1 +c2c2−d2 [d2 + (h − cos θ)2 ](2.12)

なお，この θ は原子 i を中心とした組み合わせ j − i − k の内角である．(図 2.1)．

k

θ ijk
j

i

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

実際の計算では，ある原子 i からカットオフ距離 (R + D) 内に存在する原子をリストアップし，カットオフ距離内の原子 j, k, l, … に対し，結合 i − j, i − k, i − l, … の二体間に働く力，及びそれぞれの結合から二つ，例えば結合 i − j と結合 i − k を選びその三体間によって発生する力をそれぞれ求め，これらの力のベクトルを重ね合わせる
ことにより原子 i に働く力を求めた．Tersoﬀ 型ポテンシャルでは,T2(15) および T3(16)と呼ばれる 2 種類のパラメータセットが提案されている．T2 は，表面物性を優先させてフィッティングしており， クラスターの形状を正確に表現することができるが (18) ，エネルギーが ab-initio の結果と異なることが報告されている (19) ．一方，T3 は，T2 では実験と一致しない弾性定数 C44 を修正し，結合角が 109.47◦ になるようにフィッティ
ングしているため，バルクの記述には適しているが，結合角に依存した原子挙動を示すことが報告されている.(18) 本研究では T3 を用いた.

7

Table 2.1 Parameters for silicon, gallium and arsenicSi − Si(T2)Si − Si(T3)3A [eV] 3.2647×10
1.8308×103
B [eV] 9.5373×101
4.7118×102−1
λ1 [Å ]3.23942.4799−1
λ2 [Å ]1.32581.7322−1
λ3 [Å ]1.32581.7322β3.3675×10−1 1.1000×10−6n2.2956×101
7.8734×10−1c4.83811.0039×105d2.04171.6217×101
h0.0−5.9825×10−1
R [Å ]
3.0
2.85
D [Å ]
0.2
0.15

2.4

Tersoﬀ ポテンシャルにおける力の導出

ベクトルの定義は MD のそれに従うものとする. すなわちrij ≡ ri − rjである．(rij )2 = rij · rij を微分することで，以下の関係が得られる．

rij
∂rij
=
= eij ,
∂ri
rij

∂rij
rij
=−
= −eij∂rj
rij(2.13)

原子 ij のエネルギー寄与（ボンドのエネルギー）は[

1
ϕij = fR (rij ) + (bij + bji )fA (rij )2]であるので，これの ri , rj , rk による微分

fi = −∂ϕij,∂ri

fj = −

8

∂ϕij
,
∂rj

fk = −∂ϕij∂rk(2.14)

を，j > i の全ての組について加算すれば原子に働く力が求まる．

1
fi = − [fR′ (rij ) + bij fA′ (rij )] eij − fA (rij )(2

1
fj = [fR′ (rij ) + bij fA′ (rij )] eij − fA (rij )

(

2

(

1
∂bij ∂bji
fk = − fA (rij )
+
2
∂rk
∂rk

∂bij ∂bji
+
∂ri
∂ri

∂bij ∂bji
+
∂rj
∂rj

)

(2.15))(2.16))(2.17)

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

1
∂bij
∂ζij
1
= − (1 + β n ζijn )− 2n −1 β n ζijn−1
∂rm
2
∂rm
1
∂bji
∂ζji
1
= − (1 + β n ζjin )− 2n −1 β n ζjin−1
∂rm
2
∂rm

(2.18)(2.19)

これは ri , rj , rk いずれでも同じ形になるので，−1/2fA (rij ) を含めた係数部分を
1
1
bij
1
β n ζijn−1 fA (rij )Z1 = (1 + β n ζijn )− 2n −1 β n ζijn−1 fA (rij ) =4
2 1 + β n ζijn

(2.20)と整理して ζij の微分を考える．
中心
cos θ
ik
[fc′ (rik )g(θ) ∂r+ fc (rik )g ′ (θ) ∂ ∂r
∂ζij i∑33∂ri
i=exp[λ3 (rij − rik ) ]
∂rij
∂rik
2
3
∂i
(r
−
r
)
f
(r
)g(θ)(
−
)]
+3λ
ij
ik
c
ik
3
k̸=i,j∂ri∂ri(2.21)

中心∂ζij i∑ijcos θ+ 3λ33 (rij − rik )2 fc (rik )g(θ) ∂r]=
exp[λ33 (rij − rik )3 ] [fc (rik )g ′ (θ) ∂ ∂r
∂rj
j
∂rj
k̸=i,j(2.22)cos θ
ik[fc′ (rik )g(θ) ∂r+ fc (rik )g ′ (θ) ∂ ∂r
∂ζij
3
3
∂rk
k
= exp[λ3 (rij − rik ) ]ik∂rk−3λ33 (rij − rik )2 fc (rik )g(θ) ∂r]∂rk(2.23)

9cos θ の位置ベクトルでの微分を考える．原子 j, i, k の内角であるから
rij · rik(2.24)
rij rik∂ cos θ1=(cos θeij − eik )(2.25)∂rj
rij
1∂ cos θ=(−eij + cos θeik )(2.26)∂rkrik
1
1∂ cos θ ∂ cos θ∂ cos θ=(− cos θeij + eik ) +(eij − cos θeik ) = −−(2.27)
∂ri
rij
rik
∂rj
∂rk
cos θ =

これは下図に示すようにボンド ij, ik に外向き垂直な方向である．

cosθeik
−eij

k

eik
θ
i

cosθeij
j

eij
−eik

Fig.2.2 Relationship between partial diﬀerenties of cos θ at each atom position r i ,r j ,r k

10

以下のように略記して整理すると

Z2 = exp[λ33 (rij − rik )3 ]fc′ (rik )g(θ)(2.28)

Z3 = exp[λ33 (rij − rik )3 ]fc (rik )g ′ (θ)(2.29)

Z4 = exp[λ33 (rij − rik )3 ]fc (rik )g(θ)3λ33 (rij − rik )2(2.30)

切断法で応力評価する際，2 体間として扱える部分 (作用・反作用) とそうでない部分(cos θ の微分) に分けて考える必要があるので，それらに分けて表示する．
切断法で垂直応力も生じるもの (２体間) は中心∂ζij i∑={ Z4 eij + (Z2 − Z4 )eik }∂rik̸=i,j(2.31)

中心∂ζij i∑{−Z4 eij }=∂rj
k̸=i,j(2.32)∂ζij= {−(Z2 − Z4 )eik }∂rk(2.33)

せん断応力しか生じないもの (cos θ 部分)
中心 {}∂ζij i∑
Z3
Z3
(−
cos
θe
+
e
)
+
(e
−
cos
θe
)
=
ij
ik
ij
ik
rij
rik
∂ri
k̸=i,j(2.34)中心 {}∂ζij i∑
Z3
(cos
θe
−
e
)
=
ij
ik
rij
∂rj
k̸=i,j(2.35)}∂ζij { Z3= rik (−eij + cos θeik )∂rk(2.36)

j を中心とする ζji の微分については，i と j を入れ替えた形になる．
Z2′ = exp[λ33 (rji − rjk )3 ]fc′ (rjk )g(θ)(2.37)

Z3′ = exp[λ33 (rji − rjk )3 ]fc (rjk )g ′ (θ)(2.38)

Z4′ = exp[λ33 (rji − rjk )3 ]fc (rjk )g(θ)3λ33 (rji − rjk )2(2.39)

中心∂ζji j∑={−Z4 eji }∂rik̸=i,j(2.40)

中心∂ζji j∑={Z4′ eji + (Z2′ − Z4′ )ejk }∂rjk̸=i,j(2.41)∂ζji= {−(Z2 − Z4′ )ejk }∂rk(2.42)11

また cos θ の微分部分
中心 {}∂ζji j∑
Z3′
(cos
θe
−
e
)
=
ji
jk
rji
∂ri
k̸=i,j(2.43)中心 {}∂ζji j∑
Z3′
Z3′
(−
cos
θe
+
e
)
+
(e
−
cos
θe
)
=
ji
jk
ji
jk
rji
rjk
∂rj
k̸=i,j(2.44)}∂ζji { Z3′= rjk (−eji + cos θejk )∂rk(2.45)

である.

12

2.5

速度スケーリング法

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

K=

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

(2.46)

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

設定温度が TC ，式 (2.47) より求めたある時刻の温度が 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.48)√

を

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

13

2.6

格子不安定解析

2.6.1

Born の安定基準

現在にいたるまで結晶の力学的安定に関して様々な研究がなされているが, その基礎は格子力学にある (21) ．理解を容易にするため原著と表記を変えて説明するが, ひず
み εij の均一変形を受ける結晶内部エネルギー U (εij ) は Taylor 展開により

U (εij ) = U (0) + U ′ (0)εij +U ′′ (0)εij εkl + · · ·2(2.49)と表すことができる．ここで，U ′ (0)，U ′′ (0) はそれぞれ εij による 1 次，2 次導関数を表す．結晶は εij = 0 においてつりあい状態にあるから，εij を変数とする 6 次元位相空間において，内部エネルギー曲面の勾配が 0，すなわち U ′ (0) = 0 でなければならな
い．ここである状態x における (無負荷とは限らない) 応力 σij (x) と弾性係数 Cijkl (x)はそれぞれ

σij (x) =1 ∂U (x)V (x) ∂εij

(2.50)

Cijkl (x) =1 ∂ 2 U (x)V (x) ∂εij ∂εkl

(2.51)と定義されることから, U ′ (0) = 0 の条件はx=0 で応力が 0 であることに対応する．

V(x)=はxにおける結晶の体積である．したがって式 (2.49) において 3 次以降の高次項を無視して変形すると

1
U (εij ) − U (0) = V (0)Cijkl (0)εij εkl
2(2.52)

が導かれる．すなわち，弾性係数行列 Cijkl (0) が正定でなければ，∆U =U (εij )-U (0) が負になり，x=0 よりも低いエネルギー状態が存在することになる.

14

Fig.2.3 Unit cell of fcc lattice

2.6.2

ブラベー格子を用いた安定解析

上述の議論は理解しやすいが，式 (2.49) は無負荷状態における展開であることに注意が必要である．すなわち，外力下で釣り合いにある結晶の安定性も，単純に Cijkl (x)の正値性で評価できることにはならない．そこで 1970 年代に行われた格子不安定解析 (22)–(24) では，理想均一結晶の変形を，図に示すよう単位格子を記述する 6 つのパラ
メータ a1 〜a6 で代表させることで，外力下の結晶の安定性を議論している．すなわち，
内部エネルギー U を変数 ai について展開すると式 (2.49) とまったく同じ形式になるが，ひずみによる展開では 1 次微分の項が (2.50) 式の応力となり，外力下では 0 にならないのに対して, 格子パラメータによる微分 U ′ (x) = ∂U (x)/∂U ai は釣り合い状態にある任意のxで 0 となる．したがって，外力下の結晶の安定性は U ′′ (x) = ∂ 2 U (x)/∂U ai ∂U ajの正値性で議論される．この基準を用いて初めて，「応力-ひずみ経路の不安定分岐点」の存在が特定され，計算機能が十分でない当時でも,[001] 単軸引っ張り下で現れる α-Feの bcc → bct 相変化などが解明されている (25) ．

2.6.3

弾性剛性係数による安定性解析

格子パラメータによる安定性解析では，原子の熱揺動を考慮することができず，有限温度下の結晶の安定性を議論することができない．そこで Wang，Yip，Li ら (26), (27)は，応力と弾性係数によって表される弾性剛性係数 (28) の正値性によって結晶の安定性を評価する手法を提案した．現論文 (27) では有限ひずみ，温度下の結晶の自由エネ
ルギーを展開して導出しているが，単純に弾性剛性係数が，非線形弾性体の応力-ひず

15み応答を表す物理量であることを考えれば, その正値性で結晶の安定性を評価しようと考えるのは，自然な成り行きとゆえる．具体的には，Voigt 対称性 (28) をもたせた次の弾性剛性係数 Bijkl の行列式の正負で安定性を評価する．

Bijkl = Cijkl + (σil δjk + σjl δik + σik δjl + σjk δil − σij δkl − σkl δij )/2(2.53)

ここで δ はクロネッカーのデルタである．
式 (2.51) の弾性係数は状態x における内部エネルギーにのみ依存するのに対し，弾性剛性係数は状態xの外力 (応力) の寄与が考慮される．応力，弾性係数は内部エネル
ギーの 1 次，2 次導関数で表されるが，Tersoﬀ 型ポテンシャルは 3 体項があるため Cijkl
は非常に複雑な形式となる．そこで本論文では原子応力 σila を切断法によって求め，系
の応力は

σil =

N
1 ∑
σila
N α

(2.54)と求める．また前述のように，応力 0 ではない状態においてはひずみに対する応力変化が弾性剛性係数であるので，各原子を中心に周囲の原子配置を垂直方向に ±0.1%，
せん断方向に ±0.05% の微小ひずみ摂動 (模式図 2.4) を与え，その時の応力勾配 ∆σij
から数値的に評価した．

16

i

rc

Fig.2.4 Schematic view of local strain perturbation for evaluation of local
elastic stiﬀness.

17

第3章
小規模セルによる構造不均一性の包括的検討
本章では，Tersoﬀ ポテンシャルによって表されるシリコンの，引張変形下の不安定挙動（応力急減）を分子動力学シミュレーションにより観察するとともに，そのときの各原子の Bijα について詳細に検討する．熱揺動や系の不均一性を考慮した包括的な議論とするため，小規模な周期セルによる単結晶バルクならびに粒界積層構造，横方向自由境界による単結晶ナノワイヤならびに粒界バンブー構造など様々な系について解析を行った．なお，本研究の主たる目的は，系の不安定挙動の開始と Bijα の関係解
明にあるので，不安定応答の結果生じる転位や相変態などの局所変形の詳細は議論し
ない．

3.1

解析条件

3.1.1

温度の効果

Si ダイヤモンド構造の結晶方位 ⟨100⟩, ⟨010⟩, ⟨001⟩ をそれぞれ x, y, z 軸とする座標系において，単位格子を 5 × 5 × 5 並べた立方体周期セルを用いて単結晶バルクでの解析を行った（模式図 3.1 左上）．原子数は 1000 である．温度は 1 K, 300 K, 600 K の 3
通りとし，Maxwell 速度分布に従う原子配置のゆらぎを乱数を用いてそれぞれ 8 セッ
ト準備した．応力が 0 となるようにセル辺長を制御しながら 10000 fs の分子動力学計算を行った後，∆εzz = 1.0 × 10−6 の微小ひずみ増分を毎ステップ与えることで z 方向への引張シミュレーションを行った．温度制御は速度スケーリング法により毎ステッ

18

z

L

y
Periodic boundary in x, y, z(bulk model)

x

L

(430)

_
(340)

(100)

(010)

z
x

y
Periodic boundary in z
(wire model)_(100)-(010) surface (430)-(340) surface Σ5 twist grain boundary

Fig.3.1 Simulation models.

プ行い，また横方向の垂直応力が 0 となるように横方向にもひずみを与えている．
なお，内部自由度の違いから，同一温度においても周期セルが大きいほど不安定開始ひずみが小さくなることを報告しており Tersoﬀ ポテンシャルでも 3 × 3 × 3 (216 原子), 10 × 10 × 10 (8000 原子)，20 × 20 × 20 (64000 原子) と周期セルの寸法を変えたシミュレーションを行って確認している．今回用いた 5 × 5 × 5 の周期セルでの不安定開始点は，大きいセルを用いた場合に十分収束しているとは言いがたいが，(i) ばらつ
き評価のために多数回計算を行う必要がある，(ii) 次のねじり粒界モデルとの整合性，の観点から採用している．周期セルサイズの違いによるばらつきについては次章で触
れる．

3.1.2

構造不均一性の効果

完全結晶からの構造不均一として，⃝
1 表面, ⃝
2 粒界, ⃝
3 表面+粒界, の 3 通りを考慮した．先の単結晶バルクに用いた立方体セルにおいて，z 方向にのみ周期境界を適用

19

し，x, y 方向を自由境界とすれば表面が (100) および (010) で構成されるナノワイヤとなる（模式図 3.1 左下）．また，⃝2 , ⃝3 の粒界モデルに対応して，z 軸を中心に Σ5 ねじ
り粒界の角度 36.87˚回転させて切り出した (430) および (3̄40) 表面からなるナノワイ
ヤについても検討した（同図中央）．セルの寸法ならびに原子数は (100)-(010) 表面のものと同じである．これらを上下に接合した原子数 2000 の直方体セル（同図右）に，
全方向周期境界条件を適用すると無限平面の Σ5 ねじり粒界が積層するバルク構造，z方向にのみ周期境界を適用するとバンブーナノワイヤとなり，これらで⃝2 および⃝
3 の効果を議論した．構造不均一の効果を抽出するため温度は 1 K とし，単結晶バルクの時と同様に初速度分布を変えた 8 セットを準備して初期構造緩和−引張シミュレーショ
ンを行った．なお，ひずみ速度は同じであるが，初期緩和計算は 20000 fs とした．

20

Si bulk 5x5x5 cubic cell (1000 atoms)

50

static analysis

Stress, σzz, GPa

40

30

T=1K

20
T=300K
100

T=600K

0.1

0.2

0.3

0.4

0.5

0.6

0.7

Strain, εzz

Fig.3.2 Stress–strain curves under tension (bulk).

3.2

解析結果

3.2.1

応力−ひずみ曲線

図 3.2 に，単結晶バルクの応力−ひずみ曲線の全体的な傾向 (∆εzz = 0.001 毎のプロット) を示す．図には静力学解析における応力−ひずみ曲線も破線で示している．なお，静力学解析のひずみ 0.598→0.599 で現れる不連続な応力急減は Tersoﬀ ポテンシャ
ルにおけるカットオフ関数に起因するものであり，分子動力学計算のそれとは異なる．

T = 1 K の場合，応力急減まで静力学解析の応力−ひずみ曲線と一致し，不安定を生じるひずみにはほとんどばらつきがない．温度が高くなるにつれ，静力学解析の応力−
ひずみ曲線から低応力側へのずれが引張後期に大きくなる．不安定を生じるひずみも

1 K の場合に比べてばらつくが，今回のシミュレーションでは T = 300 K の方がそのばらつきが大きい．なお，不安定後の応力−ひずみ応答が T = 1 K の場合とそれ以外で異なっているが，後で示すように det Bijα < 0 の原子数変化にも違いを生じており，異なるメカニズムで変形しているものと推測される．ただし，不安定後の変形については本論文の目的から外れるために割愛する．

21

3.2.2原子弾性剛性係数のばらつき

各温度それぞれ 1 セットについて，∆εzz = 0.001 毎に評価した全原子の Bijα の行列式

det Bijα の平均，ならびに ± 標準偏差の範囲（図が煩雑になるため ∆εzz = 0.005 もし
くは 0.01 毎に表示），そして det Bijα < 0 となった原子数を図 3.3 に示した．行列式の値は 0K の完全結晶における値で無次元化しており，応力変化を右軸にスケールを取って示している．ただし，応力急減時に原子配置が乱れると，det Bijα の平均は著しく変動し，標準偏差も桁違いとなってグラフが見にくくなることから，det Bijα の平均は応力急減が底打ちする点まで，標準偏差は応力ピークまでしか表示していない．T = 1 Kの場合，不安定を生じるまでの平均値の変化は静力学解析と完全に一致している．また，引張初期にわずかに認められた標準偏差のエラーバーが，引張によって小さくな
り，不安定変形直前ではほとんど det Bijα のばらつきがなくなっていることがわかる．
他のセットも不安定開始点まではまったく同じ挙動を示していた．一方，T = 300 K では det Bijα のばらつきが大きく，平均値も静力学解析からずれている．引張によって標準偏差のエラーバーが小さくなるのは T = 1 K の場合と同じであるが，応力急減を生じるひずみよりやや手前，標準偏差の下限が 0 程度になる付近からエラーバーの範囲が再び広がりかつ不規則になる．他の 7 セットを見ても同様の傾向を示しており，エラーバーの範囲が連続的に減少している範囲（ひずみ 0.3 程度）までは det Bijα の変化は同じであるが，エラーバーが乱れ始めるところからは det Bijα の平均変化の傾向は，8 ケースそれぞれでまちまちであった．T = 600 K の場合はパルス的に標準偏差が拡大している点が引張途中にところどころに認められる．また，標準偏差の幅が減少する傾向があるとはこのグラフでは言い難い．

22

3.2.3不安定開始点の特定

応力急減前後の原子配置データを ∆εzz = 0.00001 毎に取り直し，前節と同様に全ての原子の det Bijα の変化を調べた例を，T = 1 K の場合について図 5.13 に示す．ただし，やはり標準偏差のエラーバーをすべて表示するとグラフが見難くなるので，エラーバーは 0.001 毎に表示している．応力−ひずみ曲線を見ると，ひずみ 0.413 近傍からの応力急減による折れ曲がりの前に，応力上昇の頭打ち → なだらかな減少があり，これに対応して det Bijα のばらつきがわずかに増加していることがわかる．det Bijα < 0となった原子が出現したのは図中鉛直破線および矢印で示した点であり，応力−ひずみ曲線が折れ曲がりを示す点よりわずかに前にある．このときの det Bijα < 0 原子の数は 2 であるが，その後増加して det Bijα の平均，すなわち結晶全体の安定性が「負」になり，その点が応力−ひずみの折れ曲がり点に一致する．

T = 600 K の場合について，det Bijα の平均・標準偏差，ならびに det Bijα < 0 と判定された原子数の変化を図 3.5 に示す．ただし，極めてわずかな変化を議論した先のT = 1 K の図とはスケールが異なることに注意されたい．応力−ひずみ曲線を見ると，
やはり緩やかな応力減少の後に，図中⃝1 で示したように急激な折れ曲がり点が存在す
る．熱揺らぎによる応力−ひずみ曲線の「太さ」が，この点から無くなっていること
からも，ここが系の不安定開始点であると断定できる．厳密な定義からは，応力は不安定変形開始前と不安定変形終了後の釣り合い状態でそれぞれ長時間平均したものしか定義できないが，釣り合い状態にない瞬間瞬間の「応力」を，細かい時間間隔でモニタリングしたことで，応力ピーク後のこの折れ曲がりを見いだせたことを注記しておく．det Bijα のエラーバーの変化を見ると，T = 1 K のように連続的ではないものの，
応力ピーク→折れ曲がりの間に det Bijα のばらつきが拡大していることがわかる．ただ，ばらつきのオーダーが異なること，ならびに図 3.3(c) の det Bijα < 0 原子の変化からもわかる通り，T = 1 K の場合のように原子間のわずかな揺らぎだけに起因するものではない．図 3.5(b) 中矢印⃝2 , ⃝3 で示したひずみにおける，周期セル側面の原子配置を
図 3.6 に示す．丸印で囲った所では原子配置が乱れた部分が認められる．このような局所変形が応力急減前に生じ，かつ det Bijα < 0 の原子が漸増しているため，図 3.5(a) の23det Bijα の平均は大きく振動し，T = 1 K の時のように「折れ曲がり点＝ det Bijα の平均が負になる点」とは断定できない．ただし，平均が負になる回数は明らかに応力ピー
ク後に多くなっていることを，T = 300, 600 K の全てのケースで確認している．また，
図 3.5(b) のスケールで det Bijα < 0 原子の変化を見ると，応力ピーク→折れ曲がりの間に増加率が上昇しているように見えるが，図 3.3(c) のように εzz = 0 からの変化を通してみると，静力学解析の応力−ひずみ曲線からのずれが拡大し始める点から連続的に増加しており，応力ピーク点で明確に線引きできるものではない．また，折れ曲がり
点も det Bijα < 0 の原子が急増する点に確かに対応するが，やはり変化曲線は連続であるため線引きすることは難しい．
以上，単結晶バルクの全ての結果を総合すると (1) 応力ピークの後にわずかな応力減少があり，その間に det Bijα のばらつきが増加する，(2)(1) の後に系の不安定として応力急減を生じる，(3)T = 1 K の場合は，det Bijα の平均が負になった点が (2) に一致
する，などが明らかとなった．

24

0.6
20
0.4
Static analysis
0.2
0

10

Stress, σzz, GPa

30

0

α

Stress-strain curve

0.8

30

0.6
20
0.4
10

0.2
0

0

Number of detBij<0 atoms

200

0
0.1

0.2

0.3

0.4

0.5

30

0.6
Stress-strain curve

20

0.4
10

0.2
0

0

400

α

Number of detBij<0 atoms

400

α

α

400

0.8

-0.2

-0.2

200

Strain, εzz

0
0

0.1

0.2

0.3

0.4

0.5

200

0
0

0.1

0.2

Strain, εzz(a) T = 1 K

0.3

0.4

0.5

Strain, εzz(b) T = 300 K(c) T = 600 K

Fig.3.3 Change in the average, ± standard deviation of det Bijα and numberof det Bijα < 0 atoms (bulk).

T=1K

Si 5x5x5 cubic cell (1000 atoms)

0.02

α

31.2

Emergence of det Bij<0 atomαAverage of 6x6 determinant, ΣdetBij/N,αnormalized by detBij at perfect lattice

31
Stress-strain curve
30.8

0.01

30.6

30.4

0

0.41

0.412

0.414

30.2

Strain, εzz

Fig.3.4 Detail of unstable stress drop at T = 1 K.

25

Stress, σzz, GPa

040
Static analysis

Stress, σzz, GPa

α0.8

Si 5x5x5 cubic cell (1000 atoms)

1

40
Static analysis-0.2
Number of detBij<0 atoms

T=600KSi 5x5x5 cubic cell (1000 atoms)

1

Stress, σzz, GPa

40
Stress-strain curve

Average of 6x6 determinant, ΣdetBij/N,αnormalized by detBij at perfect latticeα

Average of 6x6 determinant, ΣdetBij/N,αnormalized by detBij at perfect lattice

T=300KSi 5x5x5 cubic cell (1000 atoms)

Average of 6x6 determinant, ΣdetBij/N,αnormalized by detBij at perfect lattice

T=1K
1

Si 5x5x5 cubic cell (1000 atoms)

T=600K

α

19

0.26

0.265

0.27

250
21

200
150

20

2
100

19
50
0
0.26

0.275

22

Strain, εzz

0.265

0.27

0.275

Strain, εzz(b) Change in number of det Bijα < 0
atoms(a) Average and standard deviation ofdet BijαFig.3.5 Detail of unstable stress drop at T = 600 K.

z
x

L

y

(a) εzz=0.265(b) εzz=0.273

Fig.3.6 Snapshots at Points ⃝2 and ⃝
3 in Fig.5(b).

26Stress, σzz, GPa

Number of detBij<0 atoms, Ninst

0

Stress, σzz, GPa

1

20

3

Stress-strain curve1
21

Si 5x5x5 cubic cell (1000 atoms)

300

22

Stress-strain curveαAverage of 6x6 determinant, ΣdetBij/N,αnormalized by detBij at perfect lattice

T=600K3.3

表面・粒界の解析結果および考察

3.3.1

表面の効果

単結晶ナノワイヤの引張応力−ひずみ曲線を図 3.7 にまとめて示す．図には先の単結晶バルクで示した完全結晶での静力学解析結果，ならびに分子動力学結果を 1 例だ
け併せて示した．{100} 表面，{430} 表面のナノワイヤいずれも静力学・完全結晶の応力ひずみ曲線より高い応力を示しているが，これは表面原子の割合が高いナノワイヤの効果と考えられる．{100} 表面と {430} 表面で応力急減を生じるひずみに差は生じているものの，バルクの 1 K の場合と同様に，初速度分布を変えても不安定開始ひず
みに大きな違いは見られない．
40

Si wire 5x5x5 cubic cell (1000 atoms), T=1K
wire (430)-(340) surface
static analysis

Stress, σzz, GPa

30

20
bulk
wire (100)-(010) surface

10

0

0.1

0.2

0.3

0.4

0.5

Strain, εzz

Fig.3.7 Stress–strain curves under tension (wire).

{100} 表面のナノワイヤにおける det Bijα の分布変化ならびに det Bijα < 0 原子数の変化を図 3.8 に示す．det Bijα のばらつきは，ナノワイヤでは表面原子と内部原子で異なる分布ピークを有すると予想されるので，標準偏差を用いるのは適切ではないが，これまでの統一した評価基準としてここでも採用した．引張によってエラーバーの幅が著しく減少し，結晶内部と表面原子の力学状態の差が小さくなっていることが示唆さ

27

れる．後述の粒界を有する系を含め，その必然性は現時点では証明できないが，先のT = 300 K の時と同じく，標準偏差の下限が 0 に達したときにばらつきが再び大きくなり始める．この変化は，図 3.9 に示したように {430} 表面のナノワイヤではより顕著である．図 3.9 の下の det Bijα < 0 の原子数変化を見ると，ひずみ 0.15 近傍から，引張前に存在した det Bijα < 0 原子が減少し，それが底打ちになるのが図中⃝1 で示したばらつきが最小となる点である．その後再び緩やかに増加するが，この⃝1 →⃝2 間の原子配置変化を見ても，図 3.6 で示したような結晶内部の乱れを見つけることは困難である．
唯一見つけた変化は，図 3.10 に示すように 1 つの表面の凹凸が強調されるような変形をしている点である．これらから推測される {430} ナノワイヤの引張におけるメカニ
ズムは，⃝
1 までは引張を駆動力とした系の均一化ならびに表面原子の安定化，⃝
1 以降
は表面部分の力学状態が，引張りに対して著しく変化する領域，と考えることができる．一方，図 3.8 の {100} ナノワイヤの引張初期における，det Bijα < 0 原子の複雑な変化であるが，εzz < 0.1 の振動している領域では，図 3.11(a) に示したように (010) お
よび (01̄0) 表面上で det Bijα < 0 原子の位置が変化している (赤く着色した原子)．一方，
ひずみ 0.1 近傍で det Bijα < 0 の原子が増加しているが，このとき図 3.11(b) に示すよ
うに det Bijα < 0 が存在する表面は (010)-(01̄0) から (100)-(1̄00) にスイッチし，かつ一様に分布して変化しなくなる（図 3.8）．ひずみ 0.2 近傍で det Bijα < 0 原子はわずかに減少するが，図 3.11(c) に示したように表面稜線上の det Bijα < 0 原子列が消滅してい
る．8 セット全てで同様の変化が見られたことから，(100) と (010) という同一指数面ながら，ダイヤモンド構造からの切り出し方による表面・エッジ部構造に起因した挙動と結論づけられる．

28

T=1K

Si nanowire (100)-(010) surfaces
40
Stress-strain curve
Static analysis(perfect lattice)

0.8

30

0.6
20
0.4
10

0.2
0

Stress, σzz, GPa

αAverage of 6x6 determinant, ΣdetBij/N,αnormalized by detBij at perfect lattice

1

0

400

α

Number of detBij<0 atoms

-0.2

200

0
0

0.1

0.2

0.3

0.4

0.5

Strain, εzz

Fig.3.8 Change in average, standard deviation of det Bijα and number ofdet Bijα < 0 atoms ((100)-(010) surface wire).

T=1K

Si nanowire (430)-(340) surfacesαNumber of detBij<0 atoms

0.8

40
Stress-strain curve
Static analysis(perfect lattice)

30

0.6
2

20

0.4
1
10

0.2
0

Stress, σzz, GPa

αAverage of 6x6 determinant, ΣdetBij/N,αnormalized by detBij at perfect lattice

1

0

-0.2
500
400
300
200
100
0

0.1

0.2

0.3

0.4

0.5

Strain, εzz

Fig.3.9 Change in average, standard deviation of det Bijα and number ofdet Bijα < 0 atoms ((430)-(3̄40) surface wire).

29

L

z
y

x

(a) εzz=0.27(b) εzz=0.40

Fig.3.10 Snapshots of side surface of (430)-(3̄40) nanowire at Points ⃝1 and⃝
2 in Fig.9.

L

z
x

y

(a) εzz=0.05(b) εzz=0.15(c) εzz=0.25

Fig.3.11 Snapshots of side surface of (100)-(010) nanowire. Red circles indicate
det Bijα < 0 atoms.

30

40

T=1K

Si grain boundary 5x5x10 cell (2000 atoms)
static analysis

Stress, σzz, GPa

30bulk (case1)
GB laminate

bamboo wire

20

10

0

0.1

0.2

0.3

0.4

0.5

Strain, εzz

Fig.3.12 Stress–strain curves under tension (Σ5 twist grain boundary).

3.3.2

粒界の効果

粒界を有するバルク・ナノワイヤの応力−ひずみ曲線を図 3.12 に示す．バルク粒界
積層構造の応力急減ひずみが，単結晶バルクのそれとほぼ一致しているが，セル内部の原子数が異なることを考えると現時点では偶然の可能性も否めない．前節の単結晶ナノワイヤに比べ，T = 1 K でも応力急減ひずみにはばらつきを生じ，特に表面＋粒界の不均一部分を持つ粒界バンブー構造で顕著である．
図 3.13 に，粒界を有するバルクにおける det Bijα の平均ならびに標準偏差，更にdet Bijα < 0 となった原子数の変化をそれぞれ示す．引張によって標準偏差のエラー
バーの幅はやはり小さくなり，標準偏差の下限が 0 近傍になったとき（図中⃝1 ）からばらつきが再び増加する．ただし，これまで見てきた例と異なるのは，ばらつき増加時に
1 より前，⃝
1
det Bijα < 0 原子が減少していることである．粒界バルク構造において，⃝
αと⃝2 の間，⃝
2 と⃝3 の間の原子配置を，det Bij < 0 の原子を赤く着色して図 3.14 に示す．
α1 までは粒界面の多くの原子が det Bij < 0
det Bijα < 0 原子は粒界近傍に存在するが，⃝αとなっているのに対し，⃝1 と⃝2 の間では限られた原子のみ det Bij < 0 となっている（図 3.14(a),(b)）．一方，応力−ひずみに変化が現れる⃝
2 と⃝3 の間では，粒界から粒内

31

へ det Bijα < 0 原子が発生している（図 3.14(c)）．これらのことから，粒界を有する系
では，(1) 均一変形の限界に達する, (2) 粒界の不均一性が拡大する, (3) 粒界から局所的変形が発生する, (4) 系の不安定挙動，というメカニズムが導かれる．なお，図 3.14(a)の上の粒界部分において，Coincident Site Lattice (CSL) 原子が det Bijα < 0 となっていることが確認できるが，粒界の不均一性拡大時にはそれらが連結された様相を呈している（図 3.14(b)）．
粒界バンブーナノワイヤにおける det Bijα の平均・標準偏差，det Bijα < 0 原子数の変化を図 3.15 に示す．粒界が存在するにも関わらず，εzz = 0.0 における det Bijα の平均は単結晶ナノワイヤより低い．このことは，粒界＋表面部分の det Bijα が著しく低下（負の値）していることを示唆している．やはり引張によってエラーバーの幅が減少し，標準偏差の下限が 0 近傍の値となった時にばらつきが再び増加するが，その挙動は {430}ナノワイヤに似ている．また「均一化」の限界ひずみはこれまでの系の中で最も小さαい．図 3.15 の⃝
1 ∼⃝
3 前後の原子配置変化を，det Bij の正負で着色して図 3.16 に示す．

粒界＋表面の効果が重畳するため簡単ではないが，図 3.15 において「均一変形」する
限界の点⃝1 に対応する図 3.16(a) と，⃝1 と⃝2 の間の点である図 3.16(b) を比較すると，楕円で囲ったように (430) 表面および粒界近傍の det Bijα < 0 原子が消滅している．このひずみ範囲では，{430} ナノワイヤでも det Bijα < 0 原子が減少していることが先の
図 3.9 から確認できる．一方，⃝2 に対応する図 3.16(c), ⃝2 と⃝3 の間の図 3.16(d) を見ると，「安定化」した {430} 表面に再び det Bijα < 0 の原子が多く現れる（楕円で囲った部
分）．⃝
3 の応力ピーク後の図 3.16(e) では，粒界から {100} 表面の結晶側に乱れを生じることで粒界近傍にくびれを生じ，これが系の不安定をもたらしたと考えられる．ただし，他の系に比べ，バンブーナノワイヤのみ応力低下があまり急峻でなく，系の不安定が明確でない．このことは，応力ピーク以降確かに応力が折れ曲がりを示し低下しているが，ひずみ制御下で不安定→局所変形による緩和が繰り返されている可能性を示唆している．応力ピーク後，det Bijα < 0 原子が階段状に増加していることもそれを支持している．

32

Si Σ5 twist grain boundary (bulk-laminate)
40

1
3
1

2

30

0.6
20
0.4
10

0.2
0

Stress, σzz, GPa

Stress-strain curve(right axis)

0.8

0Static analysis(perfect lattice)

-0.2
1000
800

α

Number of detBij atomsαAverage of 6x6 determinant, ΣdetBij/N,αnormalized by detBij at perfect lattice

T=1K

600
400
200
0
0

0.1

0.2

0.3

0.4

0.5

Strain, εzz

Fig.3.13 Change in average, standard deviation of det Bijα and number ofdet Bijα < 0 atoms (laminate structure of Σ5 twist grain boundary).

z
x

L

(430)

_
(340)

(100)

(010)

y

(a)εzz=0.28 (b)εzz=0.33 (c)εzz=0.36
Fig.3.14 Snapshots of periodic cell with Σ5 twist grain boundary. Red circles
indicate det Bijα < 0 atoms.

33Si Σ5 twist grain boundary (bamboo wire)α

Number of detBij<0 atoms

40

T=1K
3

Stress-strain curve

1.5

30

2
1

1

20

0.5

10

0

Stress, σzz, GPa

αAverage of 6x6 determinant, ΣdetBij/N,αnormalized by detBij at perfect lattice

20
Static analysis(perfect lattice)

-0.5
1000

-10

800
600
400
200
0
0

0.1

0.2

0.3

0.4

0.5

Strain, εzz

Fig.3.15 Change in average, standard deviation of det Bijα and number ofdet Bijα < 0 atoms (bamboo nanowire of Σ5 twist grain boundary).

(a)εzz=0.20(b)εzz=0.24(c)εzz=0.29(d)εzz=0.33(e)εzz=0.36

Fig.3.16 Snapshots of side surface of bamboo nanowire with Σ5 twist grainboundary.Red circles indicate det Bijα < 0 atoms.

34

第4章
薄板・ナノワイヤによる表面不均一性の議論
前章でも表面による構造不均一を議論したが，原子数が 1000，2000 と小さな系のナノワイヤでの検討であるため，表面と表面の重畳したエッジ部分の効果が大きいものと考えられる．本章では，セルの寸法を大きくして原子数を増やすとともに，ナノワ
イヤだけでなく薄板状モデルについても前章と同様の [001] 引張シミュレーションを行
い，表面の det Bijα について詳細に検討する．対象とする表面も (100) だけでなく (110)
表面についても検討する．

4.1

解析条件

図 4.1 に模式的に示すように，立方体セルの y ，z 方向に周期境界を適用した薄板モ
デル，ならびに前章でも解析したナノワイヤを対象とする．前章のセルを拡大し，単位格子を 15 × 15 × 15 並べたシミュレーションセル (図 4.2(a)，原子数 27000) により，(100) 表面薄板ならびに (100)-(010) ナノワイヤの検討を行うとともに，結晶方位を回転させて [110]，[1̄10]，[001] をそれぞれ x, y, z 軸とする座標系で，同程度の寸法に切りだした周期セル (図 4.2(b)，原子数 26460) により (110) 表面薄板と (110)-(1̄10) ナノ
ワイヤの解析を行った．温度 1K で 50000fs の緩和シミュレーションを行った後，[001]方向に, 毎ステップひずみ増分 ∆εzz =1.0×10−7 を与え引張シミュレーションを行った．

35

z

y
x

(b) Nanowire model(a)Thin plate model

8.16nm8.16nm

Fig.4.1 Schematic of simulation models．

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

8.16nm

6nm

8.1

nm

z[001]y[110]

8.09nm98.0x[110](a)(100) surface model(b)(110) surface modelFig.4.2 Atom configuration in the simulation cells．

36

4.2

解析結果

4.2.1

(100) 表面薄板の解析(100) 表面薄板の応力-ひずみ曲線，完全結晶での静力学解析結果，∆εzz =1.0×10−3毎の全原子の detBijα の分布変化ならびに det Bijα < 0 原子の原子数の変化を図 4.3 に示す．また図 4.4 にひずみ 0 および 0.12 における原子位置を detBijα ＜0 の原子を赤く着色して示す．前章の (100)-(010) ナノワイヤ同様，(100) 表面薄板でも表面の効果により静力学・完全結晶の応力-ひずみ曲線より高い応力を示している．引張により標準偏差の幅が応力急減近傍まで減少していき，表面部と結晶内部の力学状態の差が小さくなるのも前章で示したとおりである．一方，前章で示した前駆的なゆらぎ増大すなわ
ちエラーバーの幅の拡大はみられない．

detBijα ＜0 の変化を見ると，初期状態では detBijα ＜0 の原子がいくつか存在するが，これは図 4.4 に示すように表面に存在する．引張によって εzz =0.12 以降は 0 となり (図4.4(b))，応力急減時に detBijα ＜0 の原子が一気に急増する．
図 4.3 の応力急減前後の全原子の detBijα の変化を ∆εzz =1.0×10−4 毎に取り直して拡
大したものを図 4.5 に示す．グラフ下部にはピークひずみ近傍で detBijα が負となった
原子数を，detBijα の値でわけて示している．グラフ中の垂直方向の点線は局所構造変化が起きたひずみを示している．上図を見ると，応力急減前，εzz =0.391 付近から応力
ひずみ曲線がゆらいでいる．しかし detBijα の平均ならびに標準偏差の幅には変化が見られない．図 4.6 に応力にゆらぎが生じる前後のセル xy 断面図を示す．ゆらぎを生じる前は y 方向 (周期境界方向) の原子列はセル辺長に平行であるのに対し，ゆらぎを生じる領域ではわずかに偏向して周期セルの出入りを生じ凹凸を生じている．
図 4.5 上に赤丸で示した応力急減による折れ曲がり点について (a) ひとつ前のひずみ，(b) 折れ曲がり点，の原子配置を図 4.7 に示す．detBijα ＜0 の原子を赤く着色して示しており，また内部の様子も見るため，上図のセル xy 断面図のうち，detBijα >0 の原子を除いたものを下に示している．図からわかるように折れ曲がり点は detBijα ＜0の原子が内部に出現した点に一致する．

37

30

0.6
20

0.4

static analysis

0.2
10
0

Number of detBij<0 atoms

-0.2

0

6000

0.2

0.4

σzz [GPa]εzz=0.393(peak)

0.8

Tensile Stress

Average of 6x6 determinant, ΣdetBij/N,normalized by detBij at perfect lattice10

εzz=0.12

4000
2000
0
0

0.1

0.2

0.3

0.4

0.5

Tensile Strain εzz

Fig.4.3 Change in the average，± standard deviation of detBijα and the number of detBijα ＜0 atoms under tension((100) thin plate) ．
ness

ick
e th
plat
t

surfacezyx(a) εzz=0.00(b) εzz=0.12

Fig.4.4 Snapshots of (100) thin plate．Red circles indicate detBijα ＜0 atoms．

38

31
0.2
30.5
Stress strain curve

30

0

29.5

σzz [GPa]31.5
εzz=0.3956

Tensile Stress

32

0.4

ij

Average of 6x6 determinant, ΣdetBij/N,
Number of detBij<0 atoms normalized by detB at perfect lattice

32.5

Average of detBij
29
28.5

-0.2
2000

1000

detBij <0
detBij <-0.1
detBij <-0.2
detBij <-0.3
detBij <-1.0

0
0.39

0.392

0.394

0.396

Tensile Strain εzz

Fig.4.5 Detail of the unstable stress drop and change in the number of detBijα＜0 atoms under tension((100) thin plate) ．

thickness t

surfacesurface
εzz=0.390xzyεzz =0.3920

Fig.4.6 Cross-sectional view of atoms before and after the onset of stress fluctuation．

39

thickness

surfacesurface

εzz=0.3954

εzz=0.3955

αεzz=0.3955(detBIJ<0)

(a)

(b)

x
z

y

α

εzz=0.3954(detBIJ<0)

Fig.4.7 Cross-sectional view of atoms around unstable stress drop．Red circles indicate detBijα ＜0 atoms ((100)thin plate)．

40

4.2.2(110) 表面薄板の解析

図 4.8 に (110) 表面薄板の結果をまとめて示す．また図 4.9 に図 4.8(a)〜(c) で示した点における原子配置を，det Bijα < 0 原子の原子を赤く着色して示す．応力ひずみ曲線の静力学からのずれは前節の (100) 表面薄板より小さい．引張により標準偏差の幅はやはり減少するが，応力急減点より前に標準偏差の幅の減少の下げ止まりが起き，再び標
準偏差の幅がわずかに増加している．detBijα ＜0 の原子数変化も複雑な挙動を示している．すなわち，初期状態で表面部に存在した detBijα ＜0 の原子 (図 4.9(a)) は，εzz =0.178まで原子数は一定で εzz =0.178 から εzz =0.223 の間に減少し 0 となる (図 4.9(b))．その後先述のように標準偏差の幅が下げ止まり再び増加する点で，表面に det Bijα < 0 の原子が増加して，表面が覆われた状態となる (図 4.9(c))．
図 4.8 の応力急減前後の変化を拡大したものを図 4.10 に示し，また図 4.11 に応力急減前後の原子配置 (xy 断面) を detBijα ＜0 の原子を赤く着色して示した．前節の (100)
表面薄板同様，内部の様子を見るために detBijα >0 の原子を除いたものを下に付している．図 4.10 において，先述のように応力ピークより前に再び表面原子の detBijα が負となっているため，detBijα の標準偏差の幅は (100) 表面薄板に比べ大きい．また (100) 表
面薄板のような応力のゆらぎは生じていない．応力ひずみの折れ曲がりは εzz =0.3933
であり，その直後の εzz =0.3934 でさらに傾きが急峻になる．detBijα の平均の変化も対応しており，εzz =0.3933 → 0.3934 の間に detBijα の平均も負となっている．
図 4.11 をみると εzz =0.3898 で初めて表面第一層の原子以外に detBijα ＜0 の原子が現れ，応力ひずみ曲線が折れ曲がるひずみ εzz =0.3933 までに同様に表面近傍に点在している．その直後，表面部から内部に detBijα ＜0 の原子が集団的に発生して変形した．
なお，図 4.10 下の detBijα ＜0 の変化を詳細に見ると，この集団的な detBijα ＜0 の原子が発生した時に detBijα が大きな負の値をもつ原子が出現しており，この時に局所変形が発生している．

41

εzz=0.393(peak)static analysis

0.6
20

0.4
0.2
10
0

Number of detBij<0 atoms

-0.2

0

0.2

0.4

σzz [GPa]30

0.8

Tensile Stress

Average of 6x6 determinant, ΣdetBij/N,normalized by detBij at perfect lattice

1

0

6000
4000

ε zz=0.178(a)ε zz=0.00

2000(c)εzz=0.360(b)εzz=0.223

0
0

0.1

0.2

0.3

0.4

0.5

Tensile Strain εzz

Fig.4.8 Change in the average，± standard deviation of detBijα and the number of detBijα ＜0 atoms under tension((110) thin plate) ．
ss
kne ction
c
i
h
t
dire

surfacezyx(a) εzz=0.00(b) εzz=0.223(c) εzz=0.360

Fig.4.9 Snapshots of side (110) thin plate．Red circles indicate detBijα ＜0 atoms．

42

Number of detBij<0 atoms

εzz=0.3933

30.5
εzz=0.3934

0.2
Stress strain curve
Average of detBij

30

0

29.5

σzz [GPa]

0.4

Tensile Stress

Average of 6x6 determinant, ΣdetBij/N,normalized by detBij at perfect lattice

31

29

-0.2
2000

1000
detBij <0
detBij <-0.1
detBij <-0.2
detBij <-0.3
detBij <-1.0

0
0.392

0.394

Tensile Strain εzz

Fig.4.10 Detail of the unstable stress drop and change in the number of detBijα＜0 atoms under tension((110) thin plate) ．

thickness direction

surfacesurface

εzz=0.3898

εzz=0.3934

εzz=0.3933

x
α

z

y

α

εzz=0.3898(detBIJ<0)αεzz=0.3933(detBIJ<0)

εzz=0.3934(detBIJ<0)

Fig.4.11 Cross-sectional view of atoms around unstable stress drop．Red circles indicate detBijα ＜0 atoms ((110)thin plate)．
43

4.2.3(100)-(010) ナノワイヤの解析

前章で行った (100)-(010) ナノワイヤについて，セルサイズを変更して解析した結果を図 4.12 に示す．図には前章の応力-ひずみ曲線を橙色で併せて示している．セルサイズを大きくしたことにより，静力学解析からのずれは前章のそれより小さくなっているがピークひずみに大きな差はない．また (100) 表面薄板に比べると表面の原子の割合の差により初期状態の標準偏差の幅が大きくなり，ピークひずみは (100) 表面薄板の 0.3955 より小さい 0.388 となっている．
図 4.13 はひずみ 0，0.1，0.2 における detBijα ＜0 原子の分布である．前章同様，引張を与えることで初期に detBijα ＜0 が存在する表面が (100)-(1̄00)(図 4.13(a)) から (010)-(01̄0) にスイッチしている (図 4.13(b)，(c)))．図 4.12 の下図の detBijα ＜0 の初期の変化はこれに対応する．その後はほぼ一定のまま，応力急減まで変化はない (図 4.12 下)．
図 4.14 に応力急減前後の変化を拡大して示す．(100) 表面薄板同様，応力ひずみ曲
線にゆらぎが見られ，また detBijα の標準偏差の幅のパルス的な増減が応力急減前に見られる．図中に点線で示した，パルス状に標準偏差の幅が大きくなった εzz =0.3912 における原子配置を図 4.15 に示した．引張後期には断面が正方形からやや菱型に変化している．エッジ部近傍の青色で示した原子は detBijα ＜-1.0 と大きな負値を示した原子である．このように大きな負値の出現により detBijα の標準偏差の幅が大きくなっていた．図 4.14 に戻ると，この点から応力が減少し始めそのまま応力急減につながってい
くが detBijα ＜0 の原子が増加するのは ∆ε=0.0004 後の εzz =0.3916 である．
図 4.16 は図 4.14 中に点線で示した 3 つの点 (0.3901，0.3912，0.3916) における原子配置である．図 4.15 で示した大きな負値が現れたエッジではなく，矢印で示すように別の対角線方向のエッジ部から局所変化を生じて応力が低下した．先の (100) 表面薄板と同様，detBijα の平均は非常に 0 に近い値まで低下しているので，わずかなゆらぎをきっかけとした系の不安定ととらえられる．

44

0.8

30

0.6
20

static analysis

0.4
0.2

10
0

Number of detBij<0 atoms

-0.2

0

0.2

0.4

σzz [GPa]εzz=0.388(peak)stress-strain curve(1000 atoms)

Tensile Stress

Average of 6x6 determinant, ΣdetBij/N,normalized by detBij at perfect lattice

1

0

6000
4000
2000

(b)εzz =0.10(a)εzz =0.00(c)εzz =0.12

0
0

0.1

0.2

0.3

0.4

0.5

Tensile Strain εzz

Fig.4.12 Change in the average，± standard deviation of detBijα and the number of detBijα ＜0 atoms under tension((100)-(010) surface nanowire)
．

z
y

x

(a) εzz=0.00(b) εzz=0.100(c) εzz=0.120

Fig.4.13 Snapshots of side (100)-(010) surface nanowire．Red circles indicate
detBijα ＜0 atoms．

45

Number of detBij<0 atoms

Stress strain curve

ε zz=0.3916

33.5
33

σzz [GPa]

0.4

0.2
32.5

Average of detBij

32
0

Tensile Stress

Average of 6x6 determinant, ΣdetBij/N,normalized by detBij at perfect lattice

34

31.5

ε zz=0.3901
ε zz=0.3912

-0.2
2000

1000

31
detBij <0
detBij <-0.1
detBij <-0.2
detBij <-0.3
detBij <-1.0

0
0.388

0.39

0.392

0.394

Tensile Strain εzz

Fig.4.14 Detail of the unstable stress drop and change in the number of detBijα＜0 atoms under tension((100)-(010) surface nanowire)．

xzx

εzz=0.3912

α

y

z

y

εzz=0.3912(detBIJ<0)

Fig.4.15 Snapshots of side and cross-sectional view of atoms εzz =0.3912．Redcircles indicate detBijα ＜0 atoms((100)-(010) surface nanowire)．

46εzz=0.3901

εzz=0.3916

εzz=0.3912

x
α

z

y

α

εzz=0.3901(detBIJ<0)αεzz=0.3912(detBIJ<0)

εzz=0.3916(detBIJ<0)

Fig.4.16 Cross-sectional view of atoms around unstable stress drop．Red circles indicate detBijα ＜0 atoms ((100)-(010) surface nanowire)．

47

4.2.4(110)-(1̄10) ナノワイヤの解析

図 4.17 に (110)-(1̄10) ナノワイヤの解析結果を示す．(100) 表面の場合と異なり，応力ピークを示すひずみおよび応力は (110) 表面薄板とあまり変わらず，むしろわず
かに高い値を示している．図 4.18 は図中に示したポイント (a)〜(c) における原子配置である．初期状態ではどちらの表面も detBijα ＜0 の原子でおおわれているが (図4.18(a))，εzz =0.189〜0.252 の間に減少し εzz =0.252 の図 (b) では y 方向の (1̄1̄0) 表面にのみ det Bijα < 0 の原子が存在している (図 4.18(b))．図 4.17 の detBijα ＜0 の原子数変化でみるとわからないが，エラーバーが再び大きくなり始める εzz =0.282 から detBijα＜0 の原子が y 方向表面から x 方向表面にスイッチしている (図 4.18(c))．
応力急減時の挙動をこれまで同様拡大して図 4.19 に示す．応力のゆらぎはなく，応力急減前に前駆的にエラーバーの幅が拡大しているのが特徴である．さらに，(100) 表面
薄板，(100)-(010) ナノワイヤでは局所変形が発生した時にしか見られなかった detBijα＜-0.1〜-0.3 の原子も増加している．ただし detBijα ＜0 の原子数自体には大きな変化はない．
図 4.20 に応力急減前後の断面図を示す．引張後期において，detBijα ＜0 の原子が存在する表面は中心が大きくへこんでいる．模式図 4.21 に，fcc 構造とダイヤモンド構
造の単位格子を上に，ダイヤモンド構造の原子配置を z 軸方向から見たものを下に示す．下図の原子の色はそれぞれ単位格子のそれと対応している．線で囲った正方形は

x 軸，y 軸を [110]，[1̄10] にとったときのエッジ最小周期単位であり，ワイヤの各部の原子はこの最小単位のそれと一致する．実線で囲った正方形中において，黄色の 4 原子→緑の 4 原子→青の 4 原子→赤の原子と積層した構造になるため，これに起因した非対称な Poisson 収縮を生じたためにこのような表面形状になったものと推測される．
図 4.20 に戻ると，εzz =0.4041 はエラーバーの幅が前駆的に拡大し始めた点であり，応力急減の直前である εzz =0.4048 はエッジ部から y 方向表面に detBijα ＜0 の原子が増加している．なお，この時点では detBijα ＜-1.0 という大きな負の値を持つ原子は存在しない．その後，detBijα ＜0 の原子へこんだ表面の中央部分から結晶内部に出現し，同時に detBijα ＜-1.0 という大きな負の値を持つ原子が出現している．

48

30

0.8
0.6
static analysis

20

0.4
0.2
10
0

Number of detBij<0 atoms

-0.2

0

0.2

0.4

σzz [GPa]εzz=0.405(peak)Tensile Stress

Average of 6x6 determinant, ΣdetBij/N,normalized by detBij at perfect lattice

1

0

6000
ε zz=0.282

εzz=0.189

4000(c)ε zz=0.36

2000(a)ε zz=0.00(b)ε zz=0.252

0
0

0.1

0.2

0.3

0.4

0.5

Tensile Strain εzz

Fig.4.17 Change in the average，± standard deviation of detBijα and numberof detBijα ＜0 atoms under tension((110)-(1̄10) surface nanowire) ．

zyx(a) εzz=0.00(b) εzz=0.252(c) εzz=0.360

Fig.4.18 Snapshots of side (110)-(1̄10) surface nanowire．Red circles indicate
detBijα ＜0 atoms．

49

31.5

εzz=0.4041

0.2

31

εzz=0.4048

0

30.5

σzz [GPa]εzz=0.4049
Stress strain curve

Tensile Stress

0.4

ij

Average of 6x6 determinant, ΣdetB /N,
Number of detBij<0 atoms normalized by detB at perfect latticeij

32

Average of detBij

-0.2
2000

30

detBij <0
detBij <-0.1
detBij <-0.2
detBij <-0.3
detBij <-1.0

1000

0
0.402

0.404

0.406

0.408

Tensile Strain εzz

Fig.4.19 Detail of the unstable stress drop and change in the number of detBijα＜0 atoms under tension((110)-(1̄10) surface nanowire)．

εzz=0.4041

εzz=0.4049

εzz=0.4048

x
z

y

α
α

εzz=0.4041(detBIJ<0)αεzz=0.4048(detBIJ<0)

εzz=0.4049(detBIJ<0)

Fig.4.20 Cross-sectional view of atoms around unstable stress drop．Red circles indicate detBijα ＜0 atoms ((110)-(1̄10) surface nanowire)．

50

0.25abasic unit for
diamond structure

fcc unit celly[010]y[110]

x[110]

x[100]
Fig.4.21 Schematic of diamond structure and atom position at the edge of(110)-(1̄10)surface nanowire．

51

第5章
せん断変形下の不安定挙動の検討

3 章，4 章では引張変形下の不安定挙動（応力急減）に対する解析を行ってきたが，本章では異なる変形モードであるせん断変形の際の不安定挙動を観察し，そのときの各原子の Bijα について詳細に検討する．またこれまでの解析では各原子の Bijα の行列式である detBijα に着目して解析を行ってきたが，Bijα の固有値についても検討し，不安定挙動と Bijα の関係をより詳細に検討する．

5.1

解析条件

Si ダイヤモンド構造の結晶方位 [1̄10], [1̄1̄2], [111] をそれぞれ x，y ，z 軸とする座標系において，図 5.1 に示すような単位格子を x，y ，z 方向にそれぞれ 20 × 12 × 16 並べたシュミレーションセル (図 5.2) を用いて解析を行った．原子数は 46080 である．x，y方向周期境界，z 方向自由境界の条件下で 20000fs の緩和シミュレーションを行った後に，上下 3 原子層をつかみ部として固定しながら yz 方向にせん断ひずみを毎ステップ増加させてせん断シミュレーションを行った．せん断ひずみ増分は ∆εyz =1.0×10−6 /fsに固定し制御温度を 1K，300K に変化させ温度変化の影響を調べた．また制御温度を1K とした系ではせん断ひずみ増分を ∆εyz =1.0×10−7 /fs としてひずみ速度の効果も調べた．なお，温度制御は速度スケーリング法により毎ステップ行っている．

52

A

A

A

Diamond

C'

C'
Ҁ3/4a

A

A
C

C
C

B'

B

Ҁ3a

B
A

Ҁ6/2a

C

A'

A'

A
y

A

z
A

Ҁ2/2a

A

Fig.5.1 Unit cell for shear simulation．

15.06nm

A
x

B

B

C
B

y

B'

7.6

86

nm

7.987nmx [110]z [111]y [112]

Fig.5.2 Simulation cell for shear simulation．

53

5.2

解析結果および考察

5.2.1

温度の効果

T = 1 K，ひずみ増分 ∆εyz =1.0×10−6 /fs のシミュレーションにより得られた応力ひずみ関係と detBijα の平均などをまとめて図 5.3 に示した．応力はほぼ線形に上昇した後，∆εyz =0.254 でピークを示して急減している．detBijα < -1.0 の原子が初めて現れたのは ∆εyz =0.257 である．引張と同様に，せん断変形によっても標準偏差のエラー
バーの幅は小さくなる．標準偏差のエラーバーの幅はピーク時まで減少し続けており，
引張のときに見られたような，前駆的に detBijα のゆらぎが増大するような現象はせん
断では見られない．
図 5.3 の detBijα < 0 の原子数の変化を見ると，ピーク直後のひずみ εyz =0.255 の時に

det Bijα <0 となった原子がパルス状に急増している．図 5.4 に応力急減の前後の原子配置を，det Bijα < 0 の原子を赤く，つかみ部を橙色に着色して示す．応力ピーク時点には det Bijα < 0 となった原子は存在しないが，直後のひずみ εyz =0.255 の瞬間に系を横断するせん断帯状の detBijα < 0 の原子群が現れる．なお，この時点では detBijα < -1.0と大きな負の値を持つ原子は生じていない．その後，応力急減が一時停止し底打ちを示す εyz =0.257 において detBijα < -1.0 の原子が現れており，構造変化を生じていることが示唆される．この (b) → (c) の構造変化によって，detBijα < 0 の帯状組織の傾きが変化している．図 5.5 に模式的に示すように，εyz =0.255 のときには detBijα < 0 となった原子群が A → C’，C → B’，B → A’ 方向に階段状に連なっていたのに対し，ひずみ0.257 では [1̄1̄0] 方向に連続した原子群が detBijα < 0 となっている．これは εyz =0.256
から 0.257 にかけてすべりを発生したためにもたらされている．

T = 300 K で行ったシミュレーション結果を図 5.6 に示す．なお，下図の detBijα <0 の原子数変化は T = 1 K のグラフとスケールが異なる．T = 1 K に比べてピークひずみ，ピーク応力ともに小さくなり，εyz =0.246 で 20.8GPa となった．熱搖動のため

detBijα の標準偏差の幅は極めて大きく，これまでのような変形による均一化のような傾向は認められずパルス的に標準偏差が拡大している．T = 300 K の場合，detBijα <054の原子は無負荷の状態で約 600 原子存在している．また detBijα <-1.0 と大きな負値を示す原子も約 150 原子存在するが，応力ピーク後の εyz =0.249 まで変化はない．detBijα<0 の原子はひずみ 0.1 を超えたあたりから指数関数的に急増しているが応力ひずみ応
答には変化はみられない．図 5.7 に応力急減前後の原子配置を示す．瞬間瞬間の評価であるため，熱ゆらぎにより detBijα <0 原子分布は一見ランダムであり，T = 1 K のよ
うなせん断帯状の分布は確認できない (図 5.7(a))．その後図 5.7(b) をみると T = 1 Kのときと同様，ピークひずみ直後に楕円で囲んだ部分に detBijα <0 の帯状の原子群が確認でき，その後すべりを生じるとともに構造緩和が起こり，すべりを生じた部分に

30

3

20

σyz [GPa]εyz=0.254(peak)

10

Shear stress

Average of 6x6 determinant, ΣdetBij/N,normalized by detBij at perfect lattice[1̄1̄0] 方向に連続する．

Stress strain curve

2

1

0

Number of detBij<0 atoms

0

12000
εyz=0.254(peak)

8000

α

4000

detBij <-1.00
εyz=0.257
0

0.2

0.4

Shear Strain εyz

Fig.5.3 Change in the average，± standard deviation of det Bijα under shear．
And number of unstable atoms at T = 1 K．

55

z
x

(a) εyz=0.254(peak)y(c) εyz=0.257(b) εyz=0.255

Fig.5.4 Snapshots of atoms under simple shear． Red circles indicate detBijα< 0 atoms at T = 1 K．
Initial instablity band
A

A
C'

A

C

C
B'

A

A

C

A

A

C

B'

B

B
A'

A'

B'

B'

B

B
A'

C'

C

C

B

A

C'

B'

A'

A

A

C'

B'

B

A

C'

C'

After slip

A'

A'

z
y

A

εyz=0.00

AA

εyz=0.255

AAA

εyz=0.257

Fig.5.5 Schematic of instability band before and after slips．

56

A

30

20

σyz [GPa]εyz=0.246(peak)

10

Shear stress

Average of 6x6 determinant, ΣdetBij/N,normalized by detBij at perfect lattice3

Stress strain curve

2

1

0

Number of detBij<0 atoms

0
20000
16000
εyz=0.246(peak)

12000
8000

α

detBij <-1.0

4000
0
-4000

0

εyz=0.249
0.4

0.2

Shear Strain εyz

Fig.5.6 Change in the average，± standard deviation of det Bijα under shear．
And number of unstable atoms at T = 300 K．

zxy(a) εyz=0.246(peak)(b) εyz=0.248(c) εyz=0.249

Fig.5.7 Snapshots of atoms under simple shear． Red circles indicate detBijα< 0 atoms at T = 300 K．
57

5.2.2ひずみ速度の効果

T = 1 K でひずみ増分を ∆εyz =1.0×10−7 /fs とゆっくり変形させたシミュレーション
結果を図 5.8 に示す．∆εyz =1.0×10−6 /fs のピークひずみ εyz =0.254 に比べ，わずかに小さくなり εyz =0.247 でピークを示した．その他の detBijα の標準偏差の幅の減少や，

detBijα ＜0 の原子数変化などに大きな違いはない．
応力急減前後の変化を拡大したものを，∆εyz =1.0×10−6 /fs，∆εyz =1.0×10−7 /fs そ
れぞれについて図 5.9，図 5.10 に示した．また下図の detBijα ＜0 の原子数変化については，detBijα の値に応じて示している．変形速度が遅い図 5.10 をみると，ピークひず
み直後に detBijα ＜0 となった原子が初めて現れるが，単調に増加するのではなく発生→消滅→発生→消滅と振動しながら増加し，応力-ひずみにも同様の振動が認められる．
一方ひずみ速度 ∆εyz =1.0×10−6 /fs のグラフ図 5.9 では，ピーク後単調に応力は減少し，

detBijα ＜0 の原子も単調増加している ．どちらの変形速度でも，構造変化を示す点より前に detBijα ＜-0.1 という大きな負値を示す原子が前駆的に増加し，その後構造変化が発生して detBijα がさらに大きな負値をとる原子が発生する．ひずみ速度が小さい∆εyz =1.0×10−7 /fs の応答から，(1) 応力上昇が頭打ちしてゆらぎを生じる点=detBijα＜0 の原子が現れる点，(2) 構造変化を生じて応力がカタストロフ的に急減する点=負の detBijα がゆらぎながら増加し，最終的に大きな負値を生じる「不安定」点，と考えることができる．∆εyz =1.0×10−6 /fs のシミュレーションでは変化速度が大きいため (1)と (2) 間の変化が区別できないものと考えられる．
図 5.11 に上記 (1) と (2) の間の原子配置変化を detBijα ＜0 の原子を赤く着色して示した．ピークひずみでは detBijα ＜0 の原子は存在しないが (図 5.11(a))，直後に帯状の負の原子群が現れ，detBijα ＜0 の原子数の増加，減少に対応して不安定せん断帯が発生，消滅を繰り返す (図 5.11(b)〜(d))．その後 εyz =0.2508 の図 5.11(e) において detBijα ＜-0.1 の原子が増加しているが変形はまだ生じていない．その直後のひずみ εyz =0.2509
においてすべりが発生し，detBijα ＜0 の原子群が [1̄1̄0] 方向に沿った形状となる (図5.11(f))．図 5.12 は局所変形前後での原子配置を detBijα ＜0 の原子を赤く，detBijα ＜-0.1 の原子を青く着色した拡大図である．すべり面のそれぞれの結合を白の実線，赤58

の実線で示す．シリコンには shuﬄe-set，glide-set の 2 つのすべり面があり，それぞれの面にはそれぞれ 1 本，3 本の結合がある．shuﬄe-set は結合が 1 本であるため剛性が低く，せん断変形の影響を受けやすく動きやすくなっているが，glide-set は 3 本の結合により強い剛性を示し動きにくく，変形に対しての反力が大きくなると考えられる．
図 5.12(a) では glide-set の結合に detBijα ＜-0.1 という大きな負の値を持つ原子が生じ30

3Stress strain curve20

σyz [GPa]εyz=0.247(peak)

10

Shear stress

Average of 6x6 determinant, ΣdetBij/N,normalized by detBij at perfect lattice

ている．

2

1

0

Number of detBij<0 atoms

0
12000
8000
α

εyz=0.247(peak)

4000
0
-4000

detBij< -1.0

εyz=0.251
0

0.2

0.4

Shear Strain εyz

Fig.5.8 Change in the average，± standard deviation of det Bijα under shearat ∆εyz =1.0×10−7 ．

59

εyz=0.2574

Stress strain curve

23

22.5

0.2
Average of detBij

22

21.5

0

21
20.5

-0.2
Number of detBij<0 atoms

σyz [GPa]

23.5

0.4

Shear stress

Average of 6x6 determinant, ΣdetBij/N,normalized by detBij at perfect lattice

24

15000

10000

detBij <0
detBij <-0.1
detBij <-0.2
detBij <-0.3
detBij <-1.0

5000

0
0.252

0.254

0.256

0.258

Shear Strain εyz
Fig.5.9 Detail of the unstable stress drop under shear(∆εyz =1.0×10−6 /fs)．

60

Stress strain curve
εyz=0.2497

22.5
εyz=0.2487

0.2

ε yz=0.2508 ε yz=0.2509

Average of detBij

22

21.5
21

0

20.5
20

-0.2
Number of detBij<0 atoms

23

σyz [GPa]εyz=0.2492

0.4

Shear stress

Average of 6x6 determinant, ΣdetBij/N,normalized by detBij at perfect lattice

23.5

detBij <0
detBij <-0.1
detBij <-0.2
detBij <-0.3
detBij <-1.0

15000

10000

5000

0
0.246

0.248

0.25

0.252

Shear Strain εyz
Fig.5.10 Detail of the unstable stress drop under shear(∆εyz =1.0×10−7 /fs)．

61

(a) εyz=0.2481(peak)(b) εyz=0.2487
x [110](c) εyz=0.2492

z [111]y [112](d) εyz=0.2497(e) εyz=0.2508

Fig.5.11 Snapshots of atoms aroundshear(∆εyz =1.0×10−7 /fs)．

(f) εyz=0.2509

unstablestress

drop

under

shuffle-set
Glide-set
shuffle-set
Glide-set
shuffle-set
Glide-set
shuffle-set
Glide-set
shuffle-set
Glide-set(b) εyz=0.2509(a) εyz=0.2508

Fig.5.12 Magnified views of atoms around catastrophic unstable stress drop
under shear(∆εyz =1.0×10−7 /fs)．

62

5.3

固有値解析による不安定挙動の検討

これまでの解析では各原子の Bijα の行列式である detBijα に着目して不安定挙動との関係を検討してきたが，さらに詳細に系の不安定条件を明らかにするため，3 章で解析を行った T = 1 K の単結晶バルクの解析結果，せん断シミュレーションの解析結果
において 6 × 6 行列の Bijα の固有値解析を行った．
図 5.13 は 3 章で解析を行った T = 1 K の単結晶バルクの応力急減前後の挙動をさらに拡大したものである．また図 5.14 は図 5.13 で示したポイントにおける Bijα の最小の固有値 (以下 λ1 とする) が負になった原子をカウントし，横軸に固有ベクトルの成分の値をとって示したものである．detBijα ＜0 の原子が初めて現れる εzz =0.41317 では εxx ，

εyy ，εzz ，εyz に対応する変形モードに非零の成分が発生しているがその数は 0 を境に対称となっている (図 5.14(a))．しかしその直後のひずみ εzz =0.41318 では detBijα ＜0 の原子の原子数はかわらず detBijα の標準偏差の幅にも大きな変化は見られないが，εxx ，

εyy ，εzz ，εyz に対応する変形モードの非零の成分の対称性が崩れている (図 5.14(b))．
その後 detBijα ＜0 の原子の原子数は増加していくが，固有値を見ると図 5.14(c)εzx ，εxyに対応する変形モードの成分は対称性を保っており，応力ひずみ曲線もなだらかな減少を続けている．その後 εzz =0.41328 において εzx ，εxy に対応する変形モードの成分の対称性も崩れ，その直後赤丸で示すように応力ひずみ曲線の傾きが大きくなりその点が応力−ひずみの折れ曲がり点となっている．
図 5.15 は前節の ∆εyz =1.0×10−7 /fs でのせん断シミュレーションの図 5.10 の，いくつかのポイントにおける負の固有値 λ1 の固有ベクトルの成分の値ごとの分布図である．

detBijα ＜0 となった原子の原子数が発生→消滅を繰り返しながら増加している部分のεyz =0.2487，0.2492，0.2497 は，極めてわずかな変化を議論した先のバルク引張シミュ
レーションに比べるとスケールも異なり，厳密に対称とは言い難いが各成分の正負のピーク位置は 0 を境にある程度対称である (図 5.15(a)〜(c))．その後の εyz =0.2506 においても各成分の正負ピーク位置の対称性に大きな変化はみられないが (図 5.15(d))，
その次のひずみ εyz =0.2507 において，点線で囲んだように変形方向に対応する εyz の63

ピーク形状に乱れが生じている (図 5.15(e))．そして直後の detBijα ＜-0.1 と大きな負の値を持った原子が発生した εyz =0.2508 において大きく各成分の対称性が失われる (図5.15(f))．
図 5.16 はせん断シミュレーションの応力急減前後における原子配置を detBijα ，λ1 ，λ2 の値で色付けしたものである．εyz =0.2507，0.2508 において detBijα ＜0 の原子のせん断帯が現れている部分に対応して λ1 の値も低くなっている．構造変形の生じる

εyz =0.2509 では黒丸で囲んだ部分等に構造変化による緩和により detBijα が大きな正値を示す原子が出現している．一方固有値では大きな正値を示す原子は見られず，局所構造の乱れがある場合固有値は負のままである可能性を示唆している．また λ2 にも変形発生後に負の固有値が現れる．

64

Si 5x5x5 cubic cell (1000 atoms)
εzz=0.41318

31

εzz=0.41325
εzz=0.41328

ααAverage of 6x6 determinant, ΣdetBij/N,αnormalized by detBij at perfect lattice

0.015

Emergence of det Bij<0 atom
(ε zz=0.41317)

Stress-strain curve

30.9

30.8

0.01

30.7

0.005

Stress, σzz, GPa

T=1K

30.6
0.413

0.4135Strain, εzz

Fig.5.13 Detail of unstable stress drop at T = 1 K.
10
xx
yy
zz
yz
zx
xy

Number of atoms with eigenvalue λ 1<0

Number of atoms with eigenvalue λ 1<0

10

5

0
-1

-0.5

0

0.5

1

xx
yy
zz
yz
zx
xy
5

0
-1

Eigenvector corresponding to negative eigenvalue(a) εzz=0.41317

0.51(b) εzz=0.41318
xx
yy
zz
yz
zx
xy

Number of atoms with eigenvalue λ 1<0

Number of atoms with eigenvalue λ 1<0

0

10

10

5

0
-1

-0.5

Eigenvector corresponding to negative eigenvalue

-0.5

0

0.5

1

xx
yy
zz
yz
zx
xy
5

0
-1

-0.5

0

0.5

1

Eigenvector corresponding to negative eigenvalue

Eigenvector corresponding to negative eigenvalue(c) εzz=0.41325(d) εzz=0.41328

Fig.5.14 Histogram of eigenvectors of negative eigenvalue λ1 (bulk T = 1 K)．
65

4000
xx
yy
zz
yz
zx
xy

Number of atoms with eigenvalue λ 1<0

Number of atoms with eigenvalue λ 1<0

4000

2000

2000

0
-1

-0.5

0

0.5

1

-1

-0.5

1

(b) εyz=0.2492
4000

2000

xx
yy
zz
yz
zx
xy

2000

0
-0.5

0

0.5

1

0
-1

Eigenvector corresponding to negative eigenvalue

-0.5

0

0.5

1

Eigenvector corresponding to negative eigenvalue(c) εyz=0.2497(d) εyz=0.2506
4000

4000
xx
yy
zz
yz
zx
xy

Number of atoms with eigenvalue λ 1<0

Number of atoms with eigenvalue λ 1<00.5(a) εyz=0.2487
xx
yy
zz
yz
zx
xy

xx
yy
zz
yz
zx
xy

2000

2000

0
-1

0

Eigenvector corresponding to negative eigenvalue

Number of atoms with eigenvalue λ 1<0

Number of atoms with eigenvalue λ 1<00

Eigenvector corresponding to negative eigenvalue

4000

-1

xx
yy
zz
yz
zx
xy

-0.5

0

0.5

0
-1

1

-0.5

0

0.5

Eigenvector corresponding to negative eigenvalue

Eigenvector corresponding to negative eigenvalue(e) εyz=0.2507(f) εyz=0.2508

1

Fig.5.15 Histogram of eigenvectors of negative eigenvalue λ1 (under shear T = 1 K)．
66α(a) εyz=0.2507(detB )IJα(b) εyz=0.2508(detB )IJα(c) εyz=0.2509(detB )IJ

x [110]z [111]y [112](a) εyz=0.2507(λ1)(b) εyz=0.2508(λ1)(c) εyz=0.2509(λ1)(a) εyz=0.2507(λ2)(b) εyz=0.2508(λ2)(c) εyz=0.2509(λ2)

Fig.5.16 Snapshots of atoms around unstable stress drop under shear(coloredby detBijα or λ1 or λ2 value)．
67

第6章
結論

本研究は，分子動力学法を用いてバルク，表面薄板，ナノワイヤ，粒界積層構造など様々な形状のシリコンについて変形シミュレーションを行い，系の応力の急減=不安定挙動と原子弾性剛性係数 (Bijα ) の関係について検討した．以下に，得られた結果を総括する．
第 2 章では解析手法の概要を述べた．はじめに，分子動力学法の基礎方程式を示し，原子間相互作用の評価に用いた Tersoﬀ 型ポテンシャルの概要を説明した．さらに，各原子の安定性を議論する局所格子不安定解析について概説した．
第 3 章では温度，表面，粒界の効果を包括的に検討するものとして，立方体周期セルを小規模なもの (原子数 1000) に統一して単結晶バルクならびに粒界積層構造，横方
向を自由境界とした単結晶ナノワイヤならびに粒界バンブー構造について分子動力学法による [001] 方向引張シミュレーションを行い，不安定挙動 (応力急減) 発生時の各原子の Bijα について詳細に検討した．様々な得られた結果の中で，(1) 完全結晶に近い状態（T = 1 K の単結晶バルク）では，系の不安定（応力急減）の前に応力上昇の頭打ち→なだらかな低下が存在するが，その時に det Bijα の揺らぎが増大していること，(2) 熱揺らぎや構造不均一を有する系でも，det Bijα のばらつきは引張によってはじめ減少するが，応力ピークより前のひずみで再び増大すること，(3) その点では，det Bijαの標準偏差の下限がいずれの系でもほぼ 0 となっていたこと，などが重要な発見として挙げられる．
第 4 章では表面による構造不均一について，周期セルを大きくし，(100) 表面，(110)68表面の薄板およびナノワイヤについて同様の [001] 方向引張シミュレーションにより詳細に検討した．やはり引張により標準偏差の幅が減少するが，3 章で述べたようなエラーバーの拡大を生じないケースもあった．(100) 表面ナノワイヤの結果からは，応力急減直前にエラーバーがパルス的に増大するのはエッジ部に大きな負の det Bijα が出現したためであることが示唆されている．
第 5 章では，異なる変形モードであるせん断変形シミュレーションを行い，応力急減時の各原子の Bijα について同様の検討を行った．せん断変形によっても det Bijα の標準偏差のエラーバーの幅は小さくなり，ゆらぎの増大は見られなかった．∆εyz =1.0×10−6 /fsののせん断速度では，ピーク直後のひずみにおいて，内部にせん断帯状の detBijα < 0の集団的な原子群が発生し，そこにすべりを生じて [1̄1̄0] 方向に連続した detBijα < 0 の原子群を形成するのが確認された．変形速度が遅い ∆εyz =1.0×10−7 /fs のシミュレー
ションでは，応力急減前に応力-ひずみ曲線が振動するが，これに対応して detBijα ＜0の原子が発生→消滅→発生→消滅と振動しながら増加していることがわかった．この
詳細な観察結果により，(1) 応力上昇が頭打ちしてゆらぎを生じる点=detBijα ＜0 の原子が現れる点，(2) 構造変化を生じて応力がカタストロフ的に急減する点=負の detBijαがゆらぎながら増加し，最終的に大きな負値を生じる「不安定」点であることが示唆された．さらに Bijα の固有値を調べ，変形が起きる直前に Bijα の負の最小の固有値の固有ベクトルの対称性が失われていることを示した．

69

参考文献
(1) Takahiro Numazu, Yoshitada Isono, Member IEEE, and Takeshi Tanaka, Plastic
Deformation of Nanometric Single Crystal Silicon Wire in AFM Bending Test
at Intermediate Temperatures, Journal of Microelectromechanical Systems 11(2002) 125.
(2) D. R. Clarke, M. C. Kroll, P. D. Kirchner, R. F. Cook and B. J. Hockey, Phys.
Rev. Lett. 60(1988) 2156

2159.

(3) Yuhang Jing , Qingyuan Meng and Wei Zhao, Molecular dynamics simulationsof the tensile and melting behaviours of silicon nanowires, Physica E 41 (2009)685.
(4) Yuhang Jing , Qingyuan Meng , Yufei gao, Molecular dynamics simulations of
the buckling behavior of silicon nanowires under uniaxial compression, Computational Materials Science 45 (2009) 321.
(5) S. H. Park, H. J. Kim, D. B. Lee, J. S. Lee, Y. K. Choi, O. M. Kwon, Heterogeneous crystallization of amorphous silicon expedited by external forcefields: amolecular dynamics study, Superlattices and Microstructures 35 (2004) 205.
(6) S. H. Park, J. S. Kim, J. H. Park, J. S. Lee, Y. K. Choi, O. M. Kwon, Molecular
dynamics study on size-dependent elastic properties of silicon nanocantilevers,
Thin Solid Films 492 (2005) 285.
(7) 泉聡志, 原祥太郎, 酒井信介, アモルファスシリコンの表面エネルギ・表面応力,
材料

52(3), 231-234, 2003-03-15(8) 香山正憲, 山本良一, 材料界面の知的材料設計-半導体, セラミックスの粒界・界
面の原子・電子構造計算生産研究 45(10), p723-730, 1993-1070

(9) Yoshitaka Umeno, Takayuki Kitamura, Ab initio simulation on ideal shear
strength of silicon, Materials Science and Engineering B 88 (2002) 79.
(10) 楠幸久, 渋谷陽二, 冨田佳宏, 第一原理計算によるシリコン単結晶の理想せん断
強度, 熱流体系および固体系のミクロシミュレーションに関する合同シンポジウ
ム・分子動力学シンポジウム講演論文集 2000(5) (2000) 53.

(11) 大穂正史，神戸大学修士論文, (2002).
(12) 屋代如月, 北村隆行, 大谷隆一, 野田茂之, 面心立方晶の転位発生の分子動力学シ
ミュレーションおよび前駆局所ひずみ集中の格子不安定条件による特性解析, 日本機械学会論文集，A 編，64-626 (1998), 2583.

(13) 屋代如月, 北村隆行, 大谷隆一, 引張軸に対して直角方向に拘束を受けるニッケル
結晶のへき開破壊に関する原子レベル解析, 日本機械学会論文集, A 編, 63-615(1997), 2380.
(14) 屋代如月, 冨田佳宏, 転位の発生および運動時の局所格子不安定性に関する分子動力学研究, 日本機械学会論文集, A 編, 67-656 (2001), 678.

(15) Tersoﬀ, J., New empirical approach for the structure and energy of covalentsystems, Phys. Rev. B, 37–12 (1987), 6991.
(16) Tersoﬀ, J., Empirical interatomic potential for silicon with improved elastic
properties, Phys. Rev. B, 38–14 (1988), 9902.
(17) Tersoﬀ, J., Modeling solid–state shemistry : Interatomic potentials for multicomponent systems, Phys. Rev. B, 39–8 (1988), 5566.
(18) Balamane, H., Halicioglu, T. and Tiller, W. A., Phys. Rev. B, 46 (1992), 2250.
(19) 丸山茂夫, 井上知洋, 熱流体系および固体系のミクロシミュレーションに関する
合同シンポジウム・第 5 回分子動力学シンポジウム, (2000), 11.

71

(20) 上田顯, コンピューターシミュレーション, (1990), 朝倉書店.
(21) Born, M. and Huang, K., Dynamical Theory of Crystal Lattices, (1954), OxfordUP.
(22) Milstein, F., Hill, R.,Theoretical properties of cubic crystals at arbitrary
pressure-II. Shear moduli, Phys. Rev. B, 3-4B(1971), 1130.
(23) Hill, R., On the Elasticity and Stability of Perfect Crystals at Finite Strain,
Mathematical Proceedings of the Cambridge Philosophical Society, 77 (1975)
225.
(24) Thompson, J. M. T. and Shorrock, P. A., Brifurcaional Instability of An Atomic
Lattice, Journal of Mechanics and Physics of Solids, 23 (1975) 21.
(25) Milstein, F. and Farber, B., Theoretical fcc rightlaw bcc Transition under [100]Tensile Loading, Physical Review Letters, 44-4 (1980), 277.
(26) Wang, J., S.Tip, S.R.Phillpot and Wolf, D., Phys.Rev.lett., 71-25, (1993), 4182.
(27) Wang, J., Li, S.Tip, A.R.Phillpot and Wolf, D., Phys. Rev., 52-17B(1995),
12627.
(28) Wallace, D. C., Thermodynamics of Crystals, (1972), Wiley, Newyork.
(29) K. Yashiro, T. Kitamura and R.Ohtani, Atomic Simulation on Deformation
and Fracture of Nano-Single Crystal of Nickel in Tension, JSME International
Journal, Series A, 40-4, (1997), 430.

72

関連発表論文・講演論文
学術講演 単結晶シリコンの局所格子不安定性解析：不安定変形前のゆらぎ増大
藤原

正大，屋代

如月

日本材料学会 第 1 回マルチスケールマテリアルモデリングシンポジウム (第 16回分子動力学シンポジウム), 大阪大学, (2011.5)

 単結晶シリコンのせん断変形シミュレーションおよび局所格子不安定解析
藤原

正大，屋代

如月

日本機械学会 第 24 回計算力学講演会 CMD2011, 岡山大学, (2011.10)

 単結晶シリコンのせん断シミュレーション及び AES 解析
藤原

正大，屋代

如月

日本機械学会 第 25 回計算力学講演会 CMD2012, 神戸ポートアイランド南地区,(2012.10)73

論文貼り付け用ページ

74

論文貼り付け用ページ75

論文貼り付け用ページ76論文貼り付け用ページ

77

論文貼り付け用ページ

78論文貼り付け用ページ

79論文貼り付け用ページ80

謝辞
本研究を遂行するにあたり，屋代如月准教授の広い御見識と暖かい御人柄に導かれたことは，著者にとって無類の幸福でありました．また，非常に至らない部分の多い
私に，時に厳しく懇切丁寧な御指導いただき，最大限の敬意と感謝の念を表すとともに心より御礼申し上げます．浅学非才の著者が本論文を完成させるまで研究全般に渡り直接御指導を頂いたのみならず，日頃の研究生活においても多くの御助言・御協力頂きました田中克志教授に深く感謝いたします．
研究室に配属後，共に研究活動に取り組んだ池宮一繁君，住谷昂大君，竹下和也君，
中田伸哉君，西村英晃君さらに睦門賢憲さん (現，株式会社 神戸製鋼所)，山根尭君，をはじめとする表面・界面工学研究室内諸先輩方と後輩に深く感謝いたします．また，同じ材料系として議論を交わし，多くの時間を共有した固体力学研究室の同輩諸子に御礼申し上げます．
最後に，大学生活において積極的に支援して頂いた家族に心から感謝の意を表します．ありがとうございました．

2013 年 2 月
藤原81

正大