応力誘起界面不安定性のフェーズフィールドシミュレーション
Phase-field simulation of stress induced surface instability○正正高木 知弘
冨田 佳宏

(神戸大・海科)(神戸大・工)正長谷部 忠司(神戸大・工)Tomohiro TAKAKI, Faculty of Maritime Sciences, Kobe University, Kobe
Tadashi HASEBE, Faculty of Engineering, Kobe University, Kobe
Yoshihiro TOMITA, Graduate School of Science and Technology, Kobe University, Kobe
Key Words : Phase-Field Method, Surface Instability, Heteroepitaxial Growth, Self-Assembly1

緒言

新しい半導体デバイス開発のため，エピタキシャル成長を用いた量子ドットの研究が活発に行われている．この製造方法としては，絶縁物マスクでパターニングを行う選択成長，エッチングを行った加工基板を用いる方法，基板とエピタキシャル層の格子不整合を利用した自己組織化成長などがある(1) ．中でも自己組織化成長は成長前に基板を加工する必要がなく，比較的大きさの揃った量子ドットを高密度に簡単に形成できるため，今日最も注目を集めている．
自己組織化成長では，界面エネルギーと格子不整合によるひずみエネルギーのバランスによりエピタキシャル層の成長モードが決定される．その具体的な形態および成長過程を把握するためには，数値シミュレーションによる評価が不可欠である．本研究では自己組織化成長の基本的な特性を評価することを目的とし，基板上の薄膜が格子不整合によるミスフィットひずみによって形態変化する過程を再現することのできるフェーズフィールドモデルを構築し，いくつかのシミュレーションを行った．

2

解析方法

基板，薄膜，気相によって構成される系の全自由エネ
ルギー F として，次式で示す Ginzburg-Landau 型の自由エネルギー汎関数を採用する．
 a22|∇φ| dV · · · · · · · · · (1)F =f (φ, εij ) +2V
ここで，φ は固相中で 1，気相中で 0 の値をとるフェーズ
フィールド，εij はひずみテンソル，a は φ の勾配修正項．
また f (φ, εij ) はバルク内の自由エネルギー密度であり，f (φ, εij ) = fd (φ) + p (φ) fe (φ, εij ) · · · · · · · · · (2)のように次式で示すダブルウェルポテンシャル
fd (φ) と
p (φ) = φ3 10 − 15φ + 6φ2 を乗じたひずみエネルギー密度 fe (φ, εij ) の和として表される．
fd (φ) = W g (φ) · · · · · · · · · · · · · · · · · · · · · · · · · · · (3)
1fe (φ, εij ) = Dijkl (φ) εij − ε0ij εkl − ε0kl (4)22W は界面のエネルギー障壁，g (φ) = φ2 (1 − φ) ， ε0ij
は格子不整合により生じる薄膜内の初期ひずみである．
0Dijkl (φ) は気相中で零，固相中で弾性率 Dijklの値を持
つ滑らかに変化する弾性係数であり，0Dijkl (φ) = ρ (φ) Dijkl· · · · · · · · · · · · · · · · · · · · · · (5)によって表される．ここで，ρ (φ) は固相と液相の密度を表し，次の式を用い値 τ を調整することにより，φ の分布と弾性率の分布を変えることができる．

12φ − 1+ 1 · · · · · · · · · · · · · · · (6)ρ (φ) =tanh
22τ以上より，次のような時間発展 Ginzburg-Landau 方程式を得ることができる．
∂φ=∇M ∇ δF∂tδφ= M ∇2 [2φ (1 − φ) (1 − 2φ) W · · · · · · · · · · · · · (7)+30g (φ) fe (φ, εij ) − a2 ∇2 φ
δF= Dijkl εij − ε0ij ,j = σij,j = 0 · · · · · (8)δεij
ここで，フェーズフィールド φ は保存量としている．M
は mobility である．また，界面の移動に対して応力場は瞬間的に追従するとして，常に応力平衡が満足されていると仮定している．
式 (7) 中の a と W は，式 (7) の一次元問題における定常状態を考えることで，次のように界面エネルギー σ と界面幅 δ に関係付けることができる．
3δσ· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · (9)b6σbW =· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · (10)δ
界面幅 δ を λ < φ < (1 − λ) の領域とすると，
a=b = 2 tanh−1 (1 − 2λ) · · · · · · · · · · · · · · · · · · · · · (11)が得られる．
式 (7) のフェーズフィールド方程式は差分法を用い，式(8) に対しては有限要素法を用いて離散化を行いている．
また，式 (7) の時間発展は前進差分を用いている．なお，シミュレーションの過程において基板の形態は変わらな
いとし，y < ts の領域において式 (7) 中のひずみエネル
ギー項 fe (φ, εij ) = 0 としている．

3

解析条件

図 1 は解析領域，初期条件，境界条件を示している．
薄膜の初期表面位置 ysurf は正弦関数

2πx π+ysurf = α sin+ (ts + tf ) · · · · · · (12)β2を用いて表している．ここで，α は振幅，β は波長，ts は基板厚さ，tf は薄膜厚さである．また，フェーズフィー

β(a)

tf

αFilmdY

εx0 = 0εx0 = -ε mis

Gasts

yinitial surface(0,0)x(b)

SubstratedX-dXε mis(c)Fig. 1 Computational model(d)ルド φ は
b1(ysurf − y) + 1 · · · · · · · · (13)φ=tanh2δ4

解析結果

図 2 は薄膜の界面形態変化と応力 σx の分布変化．図3 は自由エネルギー F ，Fi ，Fe の変化を示している．Feは全ひずみエネルギー，Fi は全自由エネルギー F からFe を差し引いた界面の存在による過剰自由エネルギーである．図 2 の (a)〜(e) は図 3 中のそれらと対応している．
形態変化およびエネルギー変化は，図 3 中の A，B，C，D に示すように，大きく 4 つの過程に分けることができる．A は初期緩和過程および準定常過程である．ここで用いた初期表面形態では，曲率の影響が大きく図 2(a) に示すように緩和過程で平坦な界面となっている．しかし(e)
500 [MPa]-3000

Fig. 2 Morphological change and stress distribution σxA[ 10-5](a)Free energy [µJ]としている．
格子不整合による薄膜内のミスフィットひずみ εmis
は，基板格子サイズ as と薄膜格子サイズ af を用いて，εmis = (af − as ) /as として表される．シミュレーション
では基板の形態は変化しないと仮定しているため，薄膜内よりも基板内でミスフィットひずみを考慮する方が便利である．そこで，基板内の初期ひずみ ε0ij = −εmis ，薄膜内の初期ひずみ ε0ij = 0 としている．
境界条件としては，解析領域左端部の x 方向変位，下端部の y 方向変位を拘束し，右端部に x 方向の −dXεmisの強制変位を与えている．ここで，dX は解析領域長である．フェーズフィールド φ に関しては，解析領域の左右端部において周期境界，上下端部においてノイマン条件を設定している．
応力場の計算は，フェーズフィールドに比べて極端に長い時間を必要とする．そこで，計算効率を高める目的
で 100 ステップに 1 回，応力場を解く方法を採用する．
この場合，安定した解が得られることを確認している．
また，応力場は平面ひずみ問題として取り扱っている．
次に示すシミュレーション結果では，以下のパラメー
タを用いている(2) ．解析領域長 dX = 5.0 µm，解析領
域幅 dY = 0.8 µm，差分格子と有限要素サイズ dx = dy= 0.01 µm，基板厚さ ts = 0.40 µm，薄膜厚さ tf = 0.12µm，界面エネルギー σ = 1.23 J/m2 ，界面幅 δ = 4 dx，mobility M = 0.1 m5 /Js，ヤング率 E = 104 GPa，ポアソン比 ν = 0.27，時間増分 dt = 0.025 fs，ミスフィッ
トひずみ εmis = 0.02，α = 0.01 µ m，β = 0.125 µ m，
τ =0.1．なお，本シミュレーションでは界面エネルギーと弾性係数は等方性であると仮定している．

BCD(b)(c)(d)FFe(e)FiFiFe

Step No.

[ 106]

Fig. 3 Variations of free energy

ながら，ある程度時間が経過すると図 2(b) に示すように薄膜表面の一部がくぼみ，薄膜が島状の形態に変化している．B の領域では急激にひずみエネルギー Fe が低下し，形態変化がひずみエネルギーを減らすように寄与していることが分かる．この時，薄膜表面積は増加し曲率も大きくなるため界面の過剰エネルギー Fi は増加する
が，Fe の低下量の方が大きいため全自由エネルギーは低下している．図 2(b)(c)(d) のように，薄膜がある程度の島に分割されると，F と Fe の低下速度は緩やかになる．
この領域 C では図 2(d) のように複雑な形態をした島が
図 2(e) のような円形に近づく変化が見られる．図 3 の C
では Fi が変化していないため，この形態変化が各島内でのひずみエネルギーを緩和するように生じていることが分かる．最終的に領域 D の定常状態に達している．

参考文献
(1)
(2)

中嶋一雄, エピタキシャル成長のフロンティア, (2002),
共立出版.
Wang, Y. U., Jin, Y. M., Khachaturyan, A. G., Acta
Mater., 52 (2004), 81 – 92.