フェーズフィールド法に基づいたステップ・ダイナミクスの数値計算実習
三浦 均∗

2015 年 9 月 10 日

概要
本テキストは，第 39 回結晶成長討論会（2015 年 9 月 24–26 日，同志社びわこリトリートセンター，滋賀県）で実施された実習「コンピュータを使って結晶成長を観る」のテキストである。著者が作成した java アプレットを用いて，ステップ・ダイナミクス理論について視覚的・半定量的に学習することを目的とする。結晶成長理論の基本（層成長，ステップ，テラス，らせん成長，二次元核形成，など）をある程度理解している大学生・大学院生以上を対象とする。数値計算の経験は問わない。参加者は各自，java アプレットを動作させる PC を持参すること。また，事前に javaのインストールを済ませておくこと。java アプレットは当日会場にて配布する。

∗ 名古屋市立大学大学院システム自然科学研究科，e-mail: miurah@nsc.nagoya-cu.ac.jp

1

目次
はじめに

1
1.1
1.2
1.3

ステップ・ダイナミクス . . . . . . . . . . . . . . . . . . . . . . .
様々ならせん成長 . . . . . . . . . . . . . . . . . . . . . . . . . . .
本実習の狙い . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
フェーズフィールド (PF) 法

22.1
2.2
2.3
2.4
2.5

ステップの記述法 . . . . . . . . . . . . . . . . . . . . . . . . . . .
結晶面の自由エネルギー . . . . . . . . . . . . . . . . . . . . . . .
フェーズフィールド方程式 . . . . . . . . . . . . . . . . . . . . . .
直線ステップの前進速度 . . . . . . . . . . . . . . . . . . . . . . .

Java アプレット . . . . . . . . . . . . . . . . . . . . . . . . . . .
実習課題

33.1
3.2
3.3
3.43.5臨界半径 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
らせん成長 (single spiral) . . . . . . . . . . . . . . . . . . . . . .
らせん成長 (co-rotating pair) . . . . . . . . . . . . . . . . . . . .
らせん成長 (opposite pair). . . . . . . . . . . . . . . . . . . . .
らせん成長 (group of dislocations) . . . . . . . . . . . . . . . . .

233345557881010101112121 はじめに
1.1 ステップ・ダイナミクス
ステップ・ダイナミクスは，結晶成長における本質的な物理過程のひとつである。
結晶表面には原子スケールの階段（ステップ）が存在し，ここに原子や分子などの成長単位が取り込まれることでステップは前進する。多数の前進ステップが次々に結晶表面を覆っていくことで，結晶面はステップ 1 段分ずつ高さを増していく。このような成長機構を層成長といい，これに伴うステップの動きをステップ・ダイナミクスという。結晶成長速度は，結晶化駆動力，結晶の転位（欠陥）密度，不純物の有無など，様々な要因によって複雑に変化する。よって，結晶成長メカニズムを解明するには，ステップ・ダイナミクスに対する各要因の影響をひとつずつ明らかにすることが重要である。

1.2 様々ならせん成長過飽和度が比較的低いときには，二次元核形成は生じない。この場合は，結晶表面に露出したらせん転位が主要なステップ供給源となる [1, 2, 3]。らせん転位の周囲に形成される様々な成長丘を図 1 に示す。ひとつのらせん転位 (single spiral) が結晶成長に寄与している場合は，図 1a に示すように，らせん転位の周囲をステップが渦状に巻き付いたらせん成長丘を形成する。らせん転位からステップが供給されながらそれぞれのステップが前進することで，あたかも渦が一定の角速度で回転しながら外側に広がっていくように見える。単位時間あたりの渦状ステップの周回数を活動度 (activity) とよび，ステップ供給能の大きさを表す [3]。ステップ間隔 λ は二次元核の臨界半径 ρc に比例し，すなわち過飽和度 σ に反比例する [3, 4]。よって，ステップ前進速度 V の過飽和度依存性が分かれば，らせん成長丘の面成長速度 R の過飽和度依存性を求めることができる。
一方，結晶表面に複数のらせん転位が露出し，かつ，らせん転位の間隔が十分に狭い場合は，異なるらせん転位から供給されたステップ間の相互作用を考慮する必要がある。図 1b は，Burger’s ベクトルが等しいらせん転位のペア (co-rotatingpair) によるらせん成長丘である。このときの面成長速度は，転位同士が十分に離れていれば single spiral の場合と等しいが，十分に近くなると single spiral の 2 倍に漸近する [3, 5]。図 1c は，Burger’s ベクトルの符号が異なるらせん転位のペア(opposite pair) によるらせん成長丘である。このときの面成長速度は，転位同士が十分に離れていればやはり single spiral の場合と等しいが，十分に近くなるとゼロとなり，結晶成長に寄与しなくなる [3, 5]。このようならせん転位の activity はゼロであり，inactive pair と呼ばれる。最後に，図 1d に，多数のらせん転位 (group ofdislocations) によるらせん成長丘を示す。Burger’s ベクトルの等しいらせん転位が3

図 1 様 々 な ら せ ん 成 長 。(a) ら せ ん 転 位 が ひ と つ の み 露 出 し て い る 場 合 。(b)Burger’s ベクトルが等しい 2 つのらせん転位が露出している場合。(c) Burger’sベクトルの符号が異なる 2 つのらせん転位が露出している場合。(d) 多数のらせん転位が露出している場合。

一直線状に等間隔で並んでおり，かつ転位間隔が十分に狭い場合は，らせん転位の数が増えるほど面成長速度が増加し，やがて一定値に収束する [3, 5]。しかし，図 1dのように，異なる Burger’s ベクトルのらせん転位がランダムに分布するような状況では，その面成長速度を推測することは困難である。

1.3 本実習の狙い
原子間力顕微鏡 (AFM) により，一般的な結晶表面では複数のらせん転位が成長に寄与していることが観察されている [6, 7, 8, 9]。しかし，光学観察のように水平方向空間分解能が低い観察手法の場合は，必ずしも成長に寄与するらせん転位の数やその Burger’s ベクトルを直接観察できない場合も多く，その決定には詳細な解析を要する [10]。自らが実験・観察している結晶の成長機構をできるだけ正確に推測するためには，図 1 に挙げた様々ならせん成長について，その面成長速度がどのような要素にどのように依存するのかを理解しておくことが重要であろう。
本実習では，フェーズフィールド (PF) 法に基づいたステップ・ダイナミクスの数4値計算プログラム（java アプレット）を用いて，らせん成長機構を視覚的・半定量的に理解することを目的とする。図 1 に挙げたような様々ならせん成長の面成長速度が，過飽和度，らせん転位の数，その分布などの要素にどのように依存するのかを調べ，その理由について考察できるようになる。

2 フェーズフィールド (PF) 法
2.1 ステップの記述法
PF 法では，図 2 に示す通り，結晶面の二次元高さ構造をフェーズ場 ϕ(r, t) で表現する。ここで，r は結晶表面上の特定の位置を示す位置ベクトル，t は時間を表す。ϕ(r, t) はテラスにおいて一定値を取る。また，高さの異なるテラスの境界では，テラス間を繋ぐように滑らかに変化する。この薄い遷移領域がステップに対応する。
つまり，PF 法では，ステップを厚さのないシャープな境界として扱うのではなく，有限の厚さを持った境界層として扱う。この手法の最大のメリットは，ステップの位置を明示的に追跡する必要がない点である。つまり，ϕ(r, t) がそのまま結晶表面上のステップ分布を表しているので，ϕ(r, t) を “適切に”時間発展させるだけで，結晶表面上のステップ・ダイナミクスを扱うことが可能である。多数のらせん転位を同時に考慮したり，ステップ同士の合体を扱うことも容易である。
問題は，如何にして ϕ(r, t) を “適切に”時間発展させるのか，である。これには，ステップ前進速度を決める物理過程を仮定し，それを再現できるように発展方程式を導出する必要がある。具体的には，ϕ(r, t) によって定まる結晶面の自由エネルギーを定義し，それが時間とともに必ず減少するように ϕ(r, t) の時間発展の向きを与えれば良い。ステップ前進速度が溶液相からの成長ユニット直接取り込みによって決まる場合については，Miura & Kobayashi [5] によって定式化された。しかし，

Miura & Kobayashi のモデルは，ステップ前進速度を正確に再現するために大量の計算メッシュを用いるため，計算量が膨大となり，本プログラムのように数値計算結果をリアルタイムで可視化する用途には向かない。
本 java アプレットでは，Kobayashi [11] によって定式化されたモデルを採用した。以下の 2.2, 2.3 節では，ϕ(r, t) の時間発展方程式（フェーズフィールド方程式）の導出法を述べる。

2.2 結晶面の自由エネルギー
結晶面上のある領域 Ω の自由エネルギー E を，ϕ(r, t) の汎関数として，以下のよ∫ [ 2うに定義する。

E=
Ω]

ϵ
|∇ϕ|2 + F (ϕ; r) dr25(1)図 2 Schematic illustration of step description using order parameter, ϕ(r, t),where r is the position vector on the crystal surface and t is time. ϕ(r, t) takes
a constant value at each terrace and changes smoothly at steps from one side
to another within a finite width. The step dynamics on the crystal surface isdescribed by the time evolution of ϕ(r, t).

被積分関数の第 1 項目はステップの寄与，第 2 項目はステップが存在していないテラス領域の寄与を示している。関数 F (ϕ; r) は以下で与えられる。

F (ϕ; r) = F0 (ϕ − S(r))(2)

ここで，F0 (ϕ) はポテンシャル関数とよばれ，以下で与えられる。

F0 (ϕ) = − cos ϕ − µ(ϕ − sin ϕ)(3)

ここで µ は定数である。また，S(r) はシート構造関数と呼ばれ，以下で定義される。

S(r) =∑bi Arg(r − r i )(4)i

ここで bi はらせん転位 i の Burger’s ベクトルの結晶面法線方向成分，r i はらせん転位 i の結晶表面上の位置である。Arg(x) は偏角関数であり，その定義域は (−π, π]である。

6Argument function

20151050-5-10-15-201-10-0.5
x00.5-0.50.5y1 -1

図 3 (左) ポテンシャル関数 F0 (ϕ)。µ = 0.1 の場合。(右) シート構造関数 S(r)。ひ
とつのらせん転位 (bi = 1) が原点に露出している場合。

ここで，上記の式の物理的な意味を考えよう。ポテンシャル関数 F0 (ϕ) を図 3 左に示す。ϕ = 2nπ（n は整数）の位置に極小値を持つことが分かる。ここがテラスに相当する。あるテラスからひとつ上のテラスに移動すると，自由エネルギーが 2πµだけ減少する。これが結晶化駆動力に相当する。つまり，µ は結晶化駆動力の大きさを表しており，µ < 0 であれば未飽和，µ > 0 であれば過飽和であることを意味する。また，図 3 右は，シート構造関数 S(r) を図示したものである。原点にあるらせん転位を中心に一周回ると，もとの高さに戻ってこない構造を表していることが分かる。(2) 式では，S(r) を用いてテラスの基準面を連続的にズラすことで，らせん転位を導入している。このようにして得られた関数 F (ϕ; r) に，フェーズ場の勾配項，2すなわちステップの寄与 ϵ2 |∇ϕ|2 を明示的に加えたものを，ある領域 Ω にわたって積分したのが E である。

2.3 フェーズフィールド方程式
フェーズ場 ϕ(r; t) の時間発展方程式（PF 方程式）は，熱力学第二法則，すなわ
ち，「系は自由エネルギーが減少する方向に進行する」という原理に従い，自由エネ
ルギー汎関数 E から以下のように導出される [11]。

τ∂ϕ= ϵ2 ∇2 ϕ − sin(ϕ − S(r)) + µ[1 − cos(ϕ − S(r))]∂t(5)*1
ここで，τ, ϵ は定数である。

(5) 式の微分方程式を，1 辺 ∆x の正方格子上で差分化して数値積分を行なう。
メッシュ数は 100 × 100 とした。時間微分には一次精度の風上差分法を，空間微分には二次精度の中央差分法を適用した。時間刻み幅を ∆t とする。今回の計算に用いるパラメータを表 1 にまとめた。
*1 結晶表面に吸着した不純物によるステップ・ピニングの扱い方については，文献 [12] を参照のこと。

7表 1 固定計算パラメータ。

パラメータ表記

数値

時間刻み幅

∆t0.2

格子幅

∆x0.01mobility の逆数

τ

1.0定数

ϵ0.01

2.4 直線ステップの前進速度(5) 式を用いて，様々な結晶化駆動力 µ の値に対して直線ステップの前進速度 V∞を数値計算した結果を図 4 に示す。µ ≃ 0.03 − 1.0 の範囲において V∞ が µ に比例しており，近似的にV∞ = 7.5 × 10−3 µ(0.03 <∼µ<∼ 1.0)(6)で表わされることが分かる。これは，ステップ前進速度が溶液相からの成長ユニット直接取り込みによって決まる場合 [13] の依存性を表している。一方，µ <∼ 0.03 もしくは µ >
∼ 1.0 では，V∞ は (6) 式の線形関係からずれている。これは，計算に用いたメッシュ数が少なく，格子幅 ∆x を十分に小さく設定できないという計算手法上の理由に依る。以上の理由より，本 java アプレットは µ ≃ 0.03 − 1.0 の範囲での利用を推奨する。

2.5 Java アプレット
図 5 に，java アプレットによる数値計算のスナップショットを示す。画面上部のリストから，結晶化駆動力，らせん転位の数（反時計回り，時計回りそれぞれについて独立に選択可能），二次元核形成率，不純物密度を選択できる。結晶化駆動力と二次元核形成率については，計算途中でも変更が可能である。ステップが計算領域境界に到達したら，そのまま外部に流出する。計算領域外からの新たなステップ流入は考慮しない。画面上のボタンの役割や数値の意味は表 2 を参照。

8

図4

直線ステップの前進速度 V∞ の結晶化駆動力 µ 依存性。黒塗りシンボルが javaアプレットによる計算結果を示す。破線は，計算結果をおおよそでフィッティングしたもの。

図 5 Java アプレットによるらせん成長の数値計算のスナップショット。

93 実習課題
3.1 臨界半径3.1.1 課題
今回用いたポテンシャル関数（3 式）を考慮すると，臨界半径 ρc は以下で与えられる。

ρc =κ2πµ(7)

ここで，κ はステップ単位長さあたりの自由エネルギーである。このことを用いて，κ の値を推定せよ。また，解析的には κ = 8ϵ と求められるが，この式と数値計算結果を比較せよ。

3.1.2 手順
1. らせん転位数をゼロに設定して Reset ボタンを押し，らせん転位のない結晶表面を作る。

2. 核形成率を調整して，結晶表面にひとつだけ二次元島を作る。
3. 結晶化駆動力 µ を変化させ，二次元島の半径がほとんど変わらない状態を作る。
このときの二次元島の半径が，その µ における臨界半径 ρc である。
4. (7) 式を用いて κ を求める。

3.1.3 備考
• 本計算設定において，layer 単位で表示された結晶面の平均高さは，二次元島の面積に等しい。これを利用して二次元島の半径を計算できる。

• 二次元島の成長もしくは溶解が完全に停止した状態を作ることは難しい。二次元島がわずかにでも成長していれば，その駆動力における真の臨界半径はより小さいはずであり，従ってこのとき求めた κ は真の値の上限値を与える（逆に，わずかに溶解している場合は下限値を与える）。これを利用して，κ の取り得る範囲をできるだけ正確に決定すること。

3.2 らせん成長 (single spiral)
3.2.1 課題
スパイラルステップの間隔 λ が ≃ 19ρc になることを確認せよ。また，面成長速度 R の結晶化駆動力 µ 依存性を調べ，そうなる理由を考察せよ。

10

3.2.2 手順
1. らせん転位の数（反時計回り，もしくは時計回りのいずれかひとつ）を “1”に設定して Reset ボタンを押し，結晶表面にらせん転位がひとつだけ露出した状態を作る。

2. Start (Restart) ボタンを押して計算を開始させ，経過時間とそのときの結晶面の平均高さの値を定期的に記録する。

3. 横軸に経過時間，縦軸に結晶面平均高さのグラフを描き，その傾きから面成長速度 (layer/sec) を求める。同時に，ステップ間隔を測っておいて，記録する。
4. 駆動力 µ の値を変えて手順 2–3 を繰り返し，ステップ間隔 λ と面成長速度 Rの µ 依存性を求める。

3.2.3 備考
• 計算開始直後はスパイラルステップが計算領域全体に広がっていないため，面成長速度は一定値にならない。面成長速度が一定値になったときの値を測ること。

3.3 らせん成長 (co-rotating pair)3.3.1 課題
駆動力 µ が一定のとき，面成長速度 R が転位間の距離 d に対してどのように依存するかを調べよ。また，d が一定のとき，R が µ に対してどのように依存するかを調べよ。また，そのような依存性を示す理由を考察せよ。

3.3.2 手順1. 反時計回りのらせん転位の数を “2”に，時計回りのらせん転位の数を “0”に設定して Reset ボタンを押し，結晶表面に co-rotating pair が露出した状態を作る。

2. 3.2.2 節の手順 2–3 と同様にして面成長速度 R (layer/sec) を求める。この測定を，µ を一定に保ったまま Reset ボタンを押して転位間隔 d を変えて繰り返し，R の d 依存性を調べる。
3. 手順 1 により，再びらせん転位ペアが露出した状態を作る。
4. 3.2.2 節の手順 2–3 と同様にして面成長速度 R (layer/sec) を求める。この測定を，転位間隔 d を一定に保ったまま µ の値を変えて繰り返し，R の µ 依存性を調べる。このとき，Reset ボタンを押すとらせん転位の配置が変わってしまうので，Repeat ボタンを押すように注意すること。

11

3.3.3 備考• らせん転位の位置はランダムに決まる。狙った通りの配置になるとは限らないので，根気よく Reset ボタンを繰り返し押すか，潔く諦めてその配置で測定を行なうこと。

3.4 らせん成長 (opposite pair)3.4.1 課題
3.3.1 節と同じ。

3.4.2 手順1. 反時計回りのらせん転位の数を “1”に，時計回りのらせん転位の数を “1”に設定して Reset ボタンを押し，結晶表面に opposite pair が露出した状態を作る。
2. これ以降，3.3.2 節の手順 2–4 と同じ。

3.5 らせん成長 (group of dislocations)3.5.1 課題
様々ならせん転位の配置に対して，面成長速度 R の駆動力 µ 依存性を調べよ。また，R が主にどのような要素によって決定するのか検討し，その理由について考察せよ。

3.5.2 手順
1. らせん転位の数を適当に設定して Reset ボタンを押し，結晶表面にらせん転位のグループが露出した状態を作る。

2. 3.2.2 節の手順 2–3 と同様にして面成長速度 R (layer/sec) を求める。この測定を，らせん転位の配置を保ったまま µ の値を変えて繰り返し，R の µ 依存性を調べる。このとき，Reset ボタンを押すとらせん転位の配置が変わってしまうので，Repeat ボタンを押すように注意すること。また，後の考察に用いるため，計算のスナップショット画像を保存しておくこと。

3. 手順 1–2 を繰り返し，いくつかの配置に対して R(σ) を求める。

12表2Java アプレットのボタンの役割，及び，数値の意味。

ボタン or 数値

役割 or 意味Start (Restart) ボタン

計算を開始する。Stop ボタンで一時停止した計算を再開する。

Stop ボタン

計算を一時停止する。

Repeat ボタン

今実施した計算と同じ初期状態に戻す。

Reset ボタン

初期状態に戻す。ただし，初期条件（らせん転位の位置など）は新しく再設定される。

mu チョイス

結晶化駆動力 µ を選択する。計算実施中でも変更可能。

Ns(+) チョイス

反時計回りらせん転位の個数を設定する。ら
せん転位の位置はランダムに決定される。

Ns(-) チョイス

時計回りらせん転位の個数を設定する。らせん転位の位置はランダムに決定される。

Js チョイス

二次元核形成率を決定する。核形成の位置はランダムに決定される。

N imp チョイス

吸着不純物の個数を決定する。吸着位置はランダムに決定される。

time

経過時間を表示する。

height [layer]

結晶面の平均高さを layer 単位で表示する。

ϕ の平均値を ⟨ϕ⟩ とすると，⟨ϕ⟩/(2π) の値を示す。

D max

らせん転位間距離の最大値を表示する（らせん転位の総数が 1 以下の場合は 0 と表示）。

13

参考文献
[1] F. C. Frank. The influence of dislocations on crystal growth. Discuss.
Faraday Soc., 5:48–54, 1949.
[2] W. K. Burton, N. Cabrera, and F. C. Frank. Role of dislocations in crystal
growth. Nature, 163:398–399, 1949.
[3] W. K. Burton, N. Cabrera, and F. C. Frank. The growth of crystals andthe equilibrium structure of their surfaces. Philos. Trans. R. Soc. London
Ser. A, 243(866):pp. 299–358, 1951.
[4] N. Cabrera and M. M. Levine. On the dislocation theory of evaporationof crystals. Philos. Mag., 1:450–458, 1956.
[5] Hitoshi Miura and Ryo Kobayashi. Phase-field modeling of step dynamics
on growing crystal surface: direct integration of growth units to step front.
Crystal Growth & Design, 15:2165–2175, 2015.
[6] Stephen D. Durbin, Warren E. Carlson, and Misa Saros. In situ studiesof protein crystal growth by atomic force microscopy. J. Phys. D: Appl.
Phys., 26:B128–B132, 1993.
[7] K. Maiwa, M. Plomp, W. J. P. van Enckevort, and P. Bennema. AFM
observation of barium nitrate 111 and 100 faces: Spiral growth and twodimensional nucleation growth. J. Cryst. Growth, 186:214–223, 1998.
[8] A. J. Malkin, Yu. G. Kuznetsov, and A. McPherson. In situ atomic force
microscopy studies of surface morphology, growth kinetics, defect structure and dissolution in macromolecular crystallization. J. Cryst. Growth,
196:471–488, 1999.
[9] Terry A. Land and James J. De Yoreo. The evolution of growth modes
and activity of growth sources on canavalin investigated by in situ atomic
force microscopy. J. Cryst. Growth, 208:623–637, 2000.
[10] P. G. Vekilov, Yu. G. Kuznetsov, and A. A. Chernov. Interstep interaction
in solution growth; (101) ADP face. J. Cryst. Growth, 121:643–655, 1992.
[11] Ryo Kobayashi. Selected Topics on Crystal Growth, 14th International
Summer School on Crystal Growth, pages 282–291. American Institute ofPhysics, 2010.
[12] H. Miura. Phase-field modeling of step dynamics on growing crystal surface: step pinning induced by impurities. Cryst. Growth Des., 15(8):4142–4148, 2015.
[13] A. A. Chernov. Modern Crystallography III, Crystal Growth. Springer:Berlin, 1984.

14