薄膜表面形貌的相场方法模拟

中国物理学会期刊网

0 引言

薄膜材料的表面平整度是衡量薄膜质量的一个重要的指标。而薄膜质量生长的好坏，又在一定程度上决定了表面形貌薄膜的性质及其应用范围[1]。无论是与光电设备、存储设备相关的半导体薄膜还是依赖其电学性能应用的铁电氧化物薄膜，外延生长的薄膜与衬底材料之间不可能保持晶格常数的完全一致，这使得薄膜的内部存在一定程度的错配应变，当薄膜厚度较低时，薄膜的晶格常数还能与衬底保持一致，这可以称为理想化的完全匹配(fully commensurate)的晶体薄膜，但当超过某一临界厚度时，薄膜开始出现松弛(relaxation)现象[4]，同时，由于薄膜生长过程以及表面的起伏等各类因素的作用将在薄膜的内部引入缺陷及位错，从而使得薄膜的表面粗糙度增加，降低薄膜的质量。因此，为了降低表面粗糙度，提升薄膜质量，研究薄膜的生长过程，预测薄膜的表面形貌，理解错配应变以及一些辅助缓冲层技术是如何影响薄膜的形貌是十分重要的[5]。

材料的结构组织与形貌特征的研究通常处于微观范畴，材料的微观结构能够影响材料的性能，所以材料微观结构的设计，如何去优化微观组织，解决材料中的问题，对材料性能的研究尤为重要。计算模拟是现在一个很好的方法，表现出很多优秀的特点，不仅能够再现材料结构，还能够解释材料结构与材料性能的关系，进而预测部分未知材料的结构与性能以加速新材料的研发。相场计算方法，在计算材料学领域中属于一类描述介观尺度材料结构的一类有效的计算方法[7-12]。相场方法常见于模拟枝晶生长的过程[13]，近年来也在铁性智能材料方面[15]，铁电、铁磁和多铁材料方面[13]，和半导体太阳能薄膜方面[17] 等很多材料领域有所应用。由于相场方法是基于材料热力学/动力学来展开计算，这使得该计算方法具有相当现实的物理意义。使用相场方法可以避免跟踪界面[7]，从而使得相场方法在材料结构模拟方面得到了很广泛的应用，而相场方法所涉及的自由能表达式(及部分参数)可以与具体实验结果相匹配，使得相场方法在近年来同时受到广大理论工作者与实验工作者的偏爱。目前文献中已有相当多的相场计算方法相关的综述文章[7-12]与专著[18]，相场方法也已经拓展到各种不同的材料物理性质的模拟与计算。通过引入包含错配应变效应的弹性能进入热力学自由能计算，使得相场方法在研究模拟薄膜表面形貌特征方面也有了广泛的应用。本文将把焦点定位在相场方法模拟薄膜生长与薄膜表面形貌领域，首先简单的回顾了薄膜生长的几种主要的模式，着重描述了一下临界厚度问题以及薄膜表面存在的ATG不稳定性问题。随后，概括了相场方法求解此类问题的一般性过程。然后重点介绍部分是近十余年来相场方法在计算薄膜表面问题，薄膜生长问题方面的一系列研究成果，以及处理各类问题时对相场方法进行的改进与发展，并同时讨论了相场方法对薄膜中的缺陷问题和多层薄膜问题所做的一部分工作。最后给出总结以及未来相场方法在研究此类问题的突破点与展望。

1 薄膜生长

光电薄膜器件在太阳能方面不断发展，薄膜太阳能电池作为一种光电功能薄膜，可以解决能源短缺问题且无污染易于推广[20]。薄膜太阳能电池有多种，按材料可分为硅薄膜型、化合物半导体薄膜型和有机薄膜型。化合物半导体的光电特性与材料表面的平整度有很密切的联系，例如，制备表面平整的高质量InGaAs薄膜，对于InGaAs半导体光电器件有着极其重要的意义。而InGaAs与GaAs衬底之间存在晶格失配应力，制备高质量的InGaAs薄膜一直是研究的热点与难点[21]。所以，研究薄膜的生长过程，对于理解薄膜性能，具有重要的意义。

首先我们介绍薄膜生长的几种主要过程，为了能够得到单层的薄膜结构，往往可以使用分子束外延(molecular beam epitaxy, MBE)方法[22]，通过气相技术生长外延薄膜。一般认为，Ernst Bauer在1958年首先分类总结出了3种主要的薄膜生长模式[24]，如图1所示。(a)岛状生长模式，又称为Volmer-Weber生长，这种生长模式下，首先在薄膜的衬底上原子会聚集，并逐渐形成一个小岛状的结构，薄膜原子和原子之间的作用力要比原子与衬底材料之间的作用力更强，而从热力学的角度上来说，这种生长模式提供了更多的表面能，往往适用于薄膜材料与衬底材料之间晶格常数相差较大的情况下，这是由于较大的晶格常数差距将引起巨大的错配应变并产生较大的弹性能，所以岛状生长模式虽然增加了表面能，但是却能极大幅度降低薄膜材料内部的弹性能，从而使得总自由能降低。这种生长方式在半导体材料上尤为常见，如Si/Ge系统，InAs/GaAs系统中都可以观察到此类生长。这种生长方式最大的缺点在于很容易使得薄膜表面非常不平整。(b)层状生长模式，又称为Frank-Van der Merwe生长；这种生长模式下，薄膜材料首先在衬底上生长成为薄薄的一层，当这样的一种单原子薄层生长完成之后，再在第一层的基础上生长第二层原子层，这说明了薄膜材料对衬底之间的作用力更强，明显应大于薄膜原子与原子之间的作用力，这种生长方式在氧化物薄膜生长中较为常见，这种生长方式也称为二维生长。(c)层-岛状混合生长模式，又称为Stranski-Krastanov生长[25]，是由Stranski和Krastanov在1958年发现的。在这种生长模式下，薄膜材料首先被衬底材料所吸附，从而在衬底上生长成为了薄薄的几层原子层，但当超过一定的临界厚度时，在薄膜的这几层原子层上，薄膜继续生长时，开始了岛状生长的模式，从热力学的角度来说，首先薄膜的晶格常数往往是和衬底的晶格常数相匹配的，当薄膜的厚度低于临界厚度时，整体的应变能较低，但当厚度达到临界厚度时，应变能达到一个临界值，此时，薄膜的上表面开始发生松弛现象(relaxation)，松弛将很容易引入位错，位错的引入可以降低应变能，这里应变能的降低幅度一般可以远超过薄膜进行岛状生长所引入的表面能增加，使得后一阶段的生长模式进入岛状生长模式。这类生长模式应为最常见的模式，在各类薄膜生长的过程中都有观察得到。

所以说，当薄膜厚度太高时，位错和缺陷就不可避免，这时如果考虑完美晶体的情况，就需要注意到薄膜的临界厚度，对于薄膜临界厚度的讨论，最经典的两种理论来自于Matthews-Blakeslee (MB)理论[26]以及People-Bean (PB)理论[27]。其中，根据Matthews and Blakeslee’s (M&B’s)的理论指出，外延薄膜的临界厚度为

(1a)

其中hc是临界厚度，ν是泊松比，f∈mf=(af-as)/as是薄膜的错配应变，其中af 和as分别是薄膜和衬底的晶格常数。b是伯克斯矢量的尺寸，α代表了位错线与伯克斯矢量的夹角，λ代表了滑移方向与薄膜平面内垂直于滑移面与界面交线的方向的夹角。假设ν=1/3，b=4Å，cosα=1/2，cosλ=1/2。公式可以简化为

(1b)

而People-Bean理论提出外延薄膜的临界厚度是

(2a)

他们假设以GexSi1-x合金为例，b=4Å，a平均值为5.54Å。可得

(2b)

最初的MB理论和PB理论大多建立在金属-半导体薄膜体系之中，为了更准确的描述氧化物薄膜中的临界厚度，Speck和Pompe改良了MB模型[28]，同样的Maree等人改良了PB模型[29]，可以得到

这里的β表示错切参数描述缺陷核的次连续能量。改良后的模型可以更好地与氧化物类型的薄膜进行比对。从而在氧化物薄膜材料方面有一定的应用[30]。

为了更直观的演示与比较两个模型的区别，根据相关的文献[32-34]，图2所示为对数坐标下MB模型和PB模型的预测值以及相关的实验观测值。图中可以很明显的看到，两个模型均有相关的实验数据可以与理论预测值对照。而这两个模型也一直存在争议。其中MB模型是最受到广泛接受的，也是能够与实验成功对照最多的模型。PB模型受到广泛关注主要是因为一些SiGe/Si系统和InGaAs/GaAs系统的实验结果可以很好的与PB模型对照。然而，需要注意的是，很多实验，尤其是早期的实验，均存在不同程度的松弛效应，而松弛效应并未考虑在模型里[35]。

同时，我们还应关注在薄膜生长过程中，临界厚度hc处薄膜表面由于应力作用产生的波动，薄膜表面的波动很容易导致局域表面的弯曲，从而导致原子在起伏的表面上堆积形成层-岛(S-K)生长模式，这种生长模式上的岛是由于弹性不稳定性导致的，从而岛的晶格常数和块体的一致。这种弹性不稳定性往往被称为ATG不稳定性，它首先是由Asaro, Tiller以及Grinfeld (ATG)[36]所发现的，这种一致性的岛状生长是没有缺陷的，这使得这类岛状生长在纳米电子与光电子器件当中有很强的应用潜力，但这里请注意，这种应变不稳定性会随着岛间的空间减小而降低。对于单位面积岛密度较高的薄膜来说，岛与岛之间间隔较小，随着沉积的进一步发生，岛与岛之间相连时，很容易产生缺陷或者位错环，从而导致岛上出现缺陷[38]。通过线性的不稳定性分析，可以预测临界波长以及最快的生长速度[39,40]。而表面的最终形貌，需要取决于不同波长的表面波共同作用的结果。ATG稳定性在各类系统和条件下都有所报道，例如4He的固液界面[41]，Si衬底上的SiGe薄膜[42]，聚合物基体上的聚合物薄膜[43]，金属薄膜[44]等。ATG现象的产生最基本的解释可以如图3所示[45]，如果薄膜的表面是平整的，那么表面是无法松弛的，然而，一旦表面开始振荡，波峰处的应变就可以得到松弛，波谷处的应力反而会集中，这样使得波动进一步加剧，如果还同时考虑气相沉积的过程，那么原子则更有可能沉积在波峰处，同样导致波动加剧。

2 相场计算方法

为了描述在某一特定材料衬底上外延生长的薄膜，最简单的情形是引入3种不同的序参量来表示空气层(η1)，薄膜层(η2)以及衬底层(η3)，有些计算会考虑到缓冲层(η4)。ηi分别指代的是这几相在某个位置的体积百分含量，这样对于系统中的任意一处，序参量的总和应等于1，亦即

(5)

介于两个序参量ηi和ηj的界面定义为γij。一个最基本的描述形貌的相场模型，所考虑的系统的总自由能应该包括化学能，表面能和弹性能

(6)

其中，朗道形式的化学能表达展开是

(7)

这里的ωij是双势阱的最高点，V是系统的总体积。系统的总界面能需要考虑到每两个序参量的界面的总和，这样总界面能可以写成

(8)

这里的αij是梯度能系数。这里的梯度能系数，可以通过实验过程中所得到的界面能，或者第一性原理计算得到的界面能来反向推导得到。薄膜在某一特定的应力状态下的弹性能可以使用Khachaturyan的微弹性理论计算[46]

(9)

其中Cijkl是四阶弹性张量，εij是在微弹性理论的总应变，表示的是材料的本征应变。在Khachaturyan的微弹性理论中，总应变可以拆分为各向同性部分以及各向异性部分[46]。各向同性部分反映了材料的宏观平均值，而各向异性部分反映了材料内部结构或相变，缺陷等所导致的内部应变分布不均匀现象。

对于一个守恒系统，在砷化镓衬底上的InxGa1-xAs薄膜表面形态生长的时间演化是由非线性Cahn- Hilliard方程支配

(10)

考虑薄膜的生长情况下，可以在方程的右侧增加一项[47]

(11)

其中M是动力学系数，t是时间，这里的η2表示描述薄膜的序参量，t是时间，方程右边的第二项Vd nyχ描述的是薄膜的生长，其中Vd是沉积速率，χ是一个随机数，ny是表面的法线方向，也就是说

(12)

我们的模拟中，如果我们考虑一个稳定的静态情况，那么Vd这一项应该为零。仅仅当研究薄膜沉积现象时，Vd应为一个有限值，即具体的沉积速率大小。这样，计算得到各分能量，然后求得总能量Ftotal的大小，代入相场方程迭代求解，求解方程时既可以使用差分方法，或者使用半隐性傅里叶谱方法[48]，以加快计算速度并保持相对的稳定。η3是一个固定的衬底结构，并不会随时间演化，而空气相η1可以通过公式(1)计算得到。

如果考虑到材料内部成分的问题，可以同时考虑成分场c(r,t)以及相场η(r,t)两个场，来计算得到材料的相分布以及成分分布

如果考虑存在缺陷以及缺陷移动问题，一个动态的缺陷场是必不可少的，使用Gingzburg-Landau方程(Allen-Cahn)可以很好地描述缺陷场的变化

(15)

其中L是动力学系数。对于求解带有缓冲层的薄膜结构，我们通过同时求解多序参量以为核心的Cahn-Hilliard(CH) 方程：

这里需要注意到的一点是，如果采用多序参量进行计算的话，不仅仅朗道能以及界面能的能量计算较为复杂，需要更多的相关参数，更为重要的是，计算的复杂性以及误差出错的可能性也大大增加，一种较为简化的替代方式是，一般而言，多层结构是在一层完成后再添加另一层，这样相隔一层的两层材料之间可能会有接触并考虑其界面问题，但相隔两层的材料之间接触的可能性就大大降低，这样联合求解双CH方程并不断重复，就可以得到整个系统的表面以及界面结构。当然对于一些超晶格、周期性、不考虑表面(或衬底)效应的薄膜多层结构，也可以仅使用一个序参量来区分两类不同的薄膜结构。

在模拟中，大多数材料所使用的重要参数，如晶格常数，弹性系数，错配应变等可以尽可能地使用实验工作文献的数据，部分未知材料的参数可以使用相关材料的参数估算。我们以GaAs和InAs为例，其弹性常数来源于实验[49，50]。而InxGa1-xAs的弹性系数可以使用线性插值法得到。通过选择计算过程中使用的化学能系数，梯度能系数和Cahn-Hilliard方程的动力学系数，可以调整计算的空间尺寸与时间尺寸使之与实际实验过程对照。对于薄膜体系中，一般在薄膜平面内应用周期性的边界条件。计算过程中，可以使用了两种平面波来描述薄膜的表面。一种假设是表面形态是静态的正弦波，这样薄膜厚度h(x)为

(18)

其中h0是初始平均薄膜厚度，β是振幅，k是波数，和φ是初始平面波相位。而如果采用随机表面，我们认为随机表面是n个正弦波的叠加：

(19)

其中m是波型的数量，第m列波所对应的振幅，波数，和平面波相位分别为βm，km和φm。

为了表征固态薄膜的表面平整程度，使用了均方根粗糙度(RMS)，它可以定义为

(20)

其中表示平均厚度。

3 薄膜表面的相场计算方法进展

相比理论模型来说，相场方法的优势在于可以避免追踪界面，依据自由能极小值优化获得薄膜的表面界面形貌。通过界面稳定性来与实验结果进行比对。起到解释实验机理，优化实验手段的效果。在薄膜生长模式方面，尤其是在半导体材料中，岛状生长模式或者层-岛状生长模式较为常见，具体生长情况取决于生长速度与薄膜的平面内应变大小。近年来相场模拟中两类生长均有报道，与实验结果也符合较好。Liu等人[51]首先报道了外延岛状生长的薄膜表面形貌。控制生长的层数以及生长速度可以获得不同的薄膜表面形貌。Takaki等人[47]在动力学相场计算方法，引入了生长速度项，观察到了层-岛状生长模式的半导体薄膜生长。Voigt等人[52]的研究也指出，在生长过程中，部分薄膜由于应变分布不均匀等现象导致局域应变过大，断裂形成岛状-层状的混合生长，亦即SK生长模式，Szlufarska 等人[53]在InGaAs材料的模拟也观察到了SK类型的薄膜生长，说明相场计算和已知的实验报道是保持一致的。在临界厚度的预测方面相场模拟的工作较少，近年来，为了解决薄膜松弛问题，相场方法也开展了一些工作来验证经典理论。Sheng[30]和Podmaniczky[54]分别使用了经典相场动力学方法和相场晶体方法对比了相场计算得到的临界厚度与MB理论、PB理论的对比结果。对比结果显示相场计算模拟的结果与PB理论所预测的数值一致性较高。预测薄膜的ATG不稳定性正是相场方法的优势所在。相场方法因为不预先假定界面，不追踪界面等优势，可以直接模拟薄膜的表面形貌并与实验进行对比。这里需要注意的是，ATG不稳定性是由界面能和应变能两类能量相互制衡所导致的，所以，在相场的理论模型中，必须考虑应变能的作用。薄膜表面形貌及ATG不稳定性的一些细节作者在下面的3.1节中进行了详细的描述。

3.1 薄膜形貌

Khachaturyan等人应用了其微弹性理论[55]，研究了表面的不稳定性问题。通过假设表面形貌为一个随机的叠加波，观察其表面随着时间流逝而发生的形貌变化。研究显示，在存在错配应变的情况下，一个较薄的薄膜表面会首先形成沟壑，随着时间演化使得表面呈现自发组织的纳米点结构。该结构的形成来源于自由能的极小值，提升的表面能可以有效地降低材料内部的弹性能，最终形成了纳米点结构。进而通过弹性能作用引入多重正弦波叠加的随机表面，则会最终形成了纳米岛结构。

Hu等人则研究了纳米点结构的大小以及表面能/界面能对纳米点结构的影响[56]。其研究中所示，随着衬底/薄膜之间的界面能量增大，纳米点切面与衬底平面的夹角不断增大，从而影响了纳米点的形貌与尺寸，进一步的研究显示，表面角与纳米点的尺寸是成正比的，表面角的增大意味着纳米点的尺寸也在增大。Hu等人同时研究了薄膜中存在的调幅分解问题，通过调整或控制不同的薄膜的弹性系数与剪切应变，也可以很好地控制自发分解所得到的粒子的尺寸，获得的薄膜表面的自发分解图案可以更直观体现。模拟结构结果显示，可以应用错切方式(miscut)来控制粒子的大小与形貌[57]。

平面内应变对纳米点的形貌、结构是至关重要的，绝大部分的平面内应变是xy方向等同的，当然不同衬底上不同方向上生长的薄膜，xy方向上的剪切应变或许会不一致。Ni等人通过设置不同的剪切应变，来改变纳米点的形貌，以达到其应用价值[58]。在其实验中通过模拟不同方向上(x和y)的不同错配应变(0，0.02) 以及(0.02，-0.02)所观察而得出有意思的情况。如果仅一方向上存在拉伸应变，薄膜呈现波纹形状，而在一方向拉伸，另一方向压缩，则容易出现纳米点与纳米线同时存在的情况。

Tomita等人在薄膜形貌上考虑了薄膜生长的因素，在Cahn-Hilliard方程中引入沉积速率，如公式(11)[47]。研究在某一恒定沉积速率下Vd=1×10-11 m/s 下不同初始形貌、不同界面能所导致的纳米点形貌、结构以及尺寸的变化情况。Tomita阐述了在时间演化下岛状生长过程，不仅演示出了纳米点的形貌，也同时给出了x方向上正应力场的分布，还展现出在不同薄膜-衬底的界面能情况下的效果[47]。Wang等人也采用相场模拟技术，在镓砷衬底上逐步生长出In0.3Ga0.7As薄膜的过程，如图4所示[17]。

除了了解薄膜的表面形貌特征，对于不同形貌的元素分布，也会影响相关材料器件的特性。Ni等人[59]通过同时求解方程相场方程(13)、浓度场方程(14)，可以很好地描述不同形貌的纳米点结构下的浓度变化。实验采用金字塔、拱顶、截断金字塔形纳米点的浓度场分布，用来表示Ge原子偏好聚集在量子点结构的顶部、角落、尖端等地，而靠近衬底，三相交点的连接部分浓度较低[59]。

吉野雅彦等人在衬底上加工设计了一系列的凹槽，研究在此上生长金薄膜纳米点的形貌[60]。通过两者不同的过程，可以生长得到仅在平台表面的规则分布纳米点结构以及分别在平台和凹槽都有分布的纳米点结构。这样，可以通过凹槽诱导的方式来生长排列规则整齐的纳米点材料。

而Martin Grant等人开创出相场晶体方法[61]，也包含了薄膜表面模拟的机制。他提出随着薄膜的不断生长，沟壑开始出现，如果不能控制纳米点之间的间隙，空洞有可能在薄膜内部出现，从而形成缺陷。相场晶体方法的优势在于可以很直接的看到原子在薄膜表面的累积过程，这一过程可以很清晰[61]。

近年来，Wang等人在模拟薄膜生长的过程中，同时考虑了调幅分解引起的界面不稳定性问题[62]，通过引入不均衡应变来描述不同相之间的应变差值

(21)

模拟结果显示没有错配应变，随着时间推移，在薄膜中的双相自发分解过程中存在3%错配应变时以及没有宏观错配应变但存在1%相错配应变时的情形，与没有错配应变的对比情形[62]。

Wang等人继而结合实验工作研究了一类衬底上带有贴片的自组织纳米点结构[63]。实验结果表明，贴片确实有助于自组织纳米点的形成，然而贴片中是否存在应变问题。会影响到纳米点结构的形貌生长与数量。贴片带有应变的情况下，纳米点的形貌更加饱满也会与靠贴片更近一些。其实验研究电镜图也说明了这一点。

3.2 缺陷对薄膜表面的影响

有关半导体外延薄膜中的缺陷问题，尤其是位错问题，对半导体薄膜的质量来说尤为重要，在已有的薄膜相场模拟中加入缺陷或者位错的影响，可以使相场模拟的结果更加真实[65]。在计算相结构、表面结构的基础上，再把应变本身作为序参量，通过Ginzburg-Landau方程演化，可以描述应变分布，以及应变分布随时间的演化结果。Khachaturyan等人同时求解相场(方程(13))及应变场(方程(15))来获得线缺陷随着应力松弛效应的移动轨迹，他为我们清晰明确的表明错配应变在薄膜和衬底中随着时间的流逝所经过的一系类位置变化。同时还可以描述不同情况下Frank-Read源在厚衬底、薄衬底以及多层膜之间的详细的扩散情况，这些位错源随时间演化的扩散形貌各不相同[66]。

Zhuang等人通过在相场动力学方程中植入有限元网格的方法，也可以获得很好的效果[67]。可以很好地显示位错线的分布，并同时展示薄膜材料表面的位移以及变形情况，从而可以很好地预测薄膜表面的缺陷，裂纹等力学现象。还为大家直观的体现出薄膜表面在Z方向上的放大位移的情况。

Chen等人也将位错应用于薄膜与衬底之间的界面处，继而研究薄膜中的自发分解问题。采用控制变量的方法，在保持薄膜的厚度不变，使xy方向上界面处缺陷数量改变，得出不同的自发分解的结构情况。研究结果显示，薄膜中的调幅分解，很大程度上会受到界面处的错配应变的影响。在界面处的很少的缺陷，就可以引发很明显的自发分解结构的变化。表面形貌与结构可以随着缺陷的密度，薄膜的厚度而发生变化[68]。通过调整或控制薄膜的弹性系数与剪切应变，可以很好的控制自发分解所得到的粒子的尺寸。这可以应用于错切方式(miscut)来控制粒子的大小与形貌[69]。

Grant等人使用了晶体相场方法研究了带缺陷的薄膜与不带缺陷薄膜的形貌[70]，对于与衬底之间存在一定的错配应变的情况，表面粗化以及界面能增加是不可避免的，然而他们发现，带缺陷的薄膜在一定程度上可以抑制薄膜表面沟壑的生长速度，模拟出两种分别是带有缺陷以及没有缺陷的情况，事实上，在纳米岛的表面形成缺陷后，应力集中导致弹性能增大取代了表面能的增加。这一点与生长速度以及薄膜和衬底之间错配应变的关系非常密切。

Wu等人研究了处在界面，薄膜内部，以及接近薄膜表面处的缺陷对薄膜表面形貌的影响，如图5所示。很明显的接近薄膜表面的缺陷对薄膜的影响很大。缺陷的尺寸，密度等对薄膜的表面粗糙度也有不同程度的影响，带有缺陷机制的3维结构表面形貌的计算结果与实验结果相比非常类似[55]。

3.3 多层薄膜问题

在一些实际的应用体系中，多层薄膜往往是常见的，Gururajan等人最先研究了多层薄膜中的形貌界面不稳定性问题[72]。在他们的研究中，多层薄膜形貌发生弯曲，进而发生断裂是一个常见的情况。他的实验为我们呈现出的一个关于简单的三层薄膜的例子，在这个薄膜所显示的不同的形貌，一个相对较软较薄，另一个相对较硬较厚，经过实验演化，从开始不稳定状况到薄膜发生断裂的形貌的对比过程。而其中弹性各向异性、薄膜厚度、薄膜的体积百分比以及薄膜之间的间隔，都有可能会引起薄膜表面的不稳定性，进而发生薄膜断裂问题。Zaeem等人使用了有限元方法耦合相场动力学，并引入弹性能[73]。计算了5层薄膜所能展现的不同演化结构，类型Ⅰ-Ⅳ分别显示了不同条件下的多层薄膜演化结果。薄膜有可能断裂变成颗粒形貌(类型Ⅰ)，也有可能某一条消失，另一条扩展(类型Ⅱ)，或者合并(类型Ⅲ)，或者合并并扩散(类型Ⅳ)。作者可以根据不同的条件得到演化相图。

这里请注意上述所提到的多层薄膜均是两相结构的薄膜。为了解决在半导体薄膜中使用的缓冲层技术以及多层缓冲层技术，Wu等人在原有的三个序参量的基础上加入第4个序参量以描述缓冲层结构[74]。这样通过联解双Cahn-Hilliard方程可以很好的分析带有缓冲层结构的多层薄膜体系的应力、应变、表面等等各类问题。研究显示，随着缓冲层厚度的增加，薄膜的表面粗糙度不断的降低，如图6所示[75]。

这个模型也可以进一步的拓展成多序参量结构，用来模拟量子阱、量子点结构或者并设计新一代的缓冲层结构、多层缓冲层结构以及梯度缓冲层结构，或进一步研究缓冲层的缓冲机制。

4 结语

薄膜表面以及薄膜生长决定了薄膜质量，在相场方法中考虑弹性能以及弹性效应，可以很好地描述薄膜的表面形貌，内部的应力应变分布，一定情况下还可以描述一些原子浓度，相分离，缺陷影响等特征。本文首先介绍了薄膜生长的基础，包括生长模式、临界厚度以及ATG表面不稳定性。通过基本的相场计算模型，重点介绍了近年来相场方法与经典理论的对照，及其应用在研究薄膜表面、缺陷薄膜作用以及多层薄膜情况。从已有的研究中我们可以看出，相场方法已经可以很好地与实验相比对，描述薄膜的形貌以及相关的物理性能，未来相场方法在模拟薄膜方面的发展应侧重于应用和预测材料性质方面：

4.1 多场的联合效应

从上述的研究工作我们可以看出，联合多场进行相场动力学计算，是当前相场计算发展的一大趋势。多场联合计算的好处更加全面、详细，而且考虑到不同场之间的相互作用时会更加接近真实情况。但是缺点也非常明显，其一就是计算量太大，所以多场方程联合求解如何优化计算是一个主要问题。其二是程序开发周期或许会延长，一些界面友好的相场计算软件(如μ-pro[76])，或许可以较好的解决这类问题。

4.2 基因组计算

对于材料应用来说，开发新型材料或者设计新型材料结构尤为重要，然而，完全通过实验手段开发新材料自然需要很高的成本，通过材料计算的方法进行一个大致方向上的预测，就显得尤为重要。相场动力学方法的一大优势在于可以很好的结合实验工作，获取实验参数或者根据第一性原理计算获得基本参数，从而实现与实验结果的高度匹配。同时相场方法还可以显示结构，耦合缺陷等一系列材料中的问题，并可以从统计意义上给出宏观性质。根据一系列参数进行高通量计算建立相关基因组数据库，包含材料结构，宏观性质，缺陷分布等是相场方法在近年来的首要任务。

4.3 机器学习的应用

基因组计算的目的是为了更快的找到适合的材料。而使用近年来兴起的机器学习技术，可以更快的更方便的找到新型材料。可以有两种途径达到寻找适合材料或者新材料的目的。(1)通过相场方法，建立材料基因组数据库，通过机器学习，满足适合条件的方法在数据库中找到新型材料。(2)从某一特殊点开始，通过相场方法计算出材料的性质，衡量与目标材料性质的差距，通过机器学习的方式，改变或调整初始条件，慢慢逼近目标材料所需要的特性。当然前一种方法的耗时较短，覆盖范围大，但建基因组数据库工程量大，后一种方法目的性强，灵活性强，不需要预先建立数据库。这两类方法都有自己的特点和优势，可以应用在不同场合。

参考文献

[1] JOHN V. Introduction to surface and thin film processes[M]. Cambridge: Cambridge University Press, 2000.

[2] ALBERTO P, VILLAIN J. Physics of crystal growth[M]. Cambridge: Cambridge University Press, 1998.

[3] OHRING M. Materials Science of Thin Films (2nd ed.)[M]. Boston: Academic Press, 2001.

[4] MARTIN L W, CHU Y H, RAMESH R. Advances in the growth and characterization of magnetic, ferroelectric, and multiferroic oxide thin films[J]. Materials Science and Engineering R, 2010, 68: 89.

[5] BRUNNER K. Si/Ge nanostructures[J]. Reports on Progress in Physics, 2002, 65: 27.

[6] STANGL J, HOLY V, BAUER G. Structural properties of self-organized semiconductor nanostructures[J]. Reviews of Modern Physics, 2004, 76: 725.

[7] CHEN L Q. Phase-field models for microstructure evolution[J]. Annual Review of Materials Research, 2002, 32: 113.

[8] STEINBACH I. Phase-field models in materials science[J]. Modelling and Simulation in Materials Science and Engineering, 2009, 17: 073001.

[9] SINGER-LOGINOVA I, SINGER H M. The phase field technique for modeling multiphase materials[J]. Reports on Progress in Physics, 2008, 71: 106501.

[10] WANG Y, LI J. Phase field modeling of defects and deformation[J]. Acta Materialia 2010, 58: 1212.

[11] CHEN L Q. Phase-field method of phase transitions/domain structures in ferroelectric thin films: A review[J]. Journal of the American Ceramic Society, 2008, 91(6): 1835.

[12] 陈龙庆. 相场模拟与材料基因组计划[J]. 科学通报, 2013, 58:3638-3641.

CHEN L Q. Phase-field method and materials genome initiative(MGI)[J].Chin Sci Bull. 2013, 58: 3638-3641. (in Chinese)

[13] 方鹏均，刘芯宇，张雷，等. 边界热通量作用下枝晶生长相场模型的有限元法模拟[J]．功能材料，2014, 45(18): 18042-18046．

FANG P, LIU X, ZHANG L, et al．Phase field simulation of dendritic growth with boundary heat flux by finite element method[J]．Journal of Functional Materials，2014, 45(18): 18042-18046．(in Chinese)

[14] 王智平，王宝成，肖荣振, 等. Si-Ni合金中小晶面枝晶生长的相场法模拟[J]．功能材料，2014, 45(18): 18042-18046．

WANG Z, WANG B, XIAO R, et al．Phase-field method simulation of faceted dendrite growth in Si-Ni alloy[J]．Journal of Functional Materials，2014, 45(18): 18042-18046. (in Chinese)

[15] 王栋, 王云志, 李巨. 铁性智能材料的计算机模拟进展[J]. 中国材料进展, 2012, 31(3): 8-14.

WANG D, WANG Y, LI J. Progress in Computer Simulations of Ferroic Smart Materials[J]. Materials China, 2012, 31(3): 8-14. (in Chinese)

[16] 王杰, 李欣凯, 刘畅, 等. 材料微结构演化的相场模拟[J]. 固体力学学报, 2016, 37(1): 1-33.

WANG J, LI X, LIU C, et al. Phase Field Simulations of Microstructure Evolution[J]. Chinese Journal of Solid Mechanics, 2016, 37(1): 1-33. (in Chinese)

[17] 王冠，翁燕华，吴平平. 相场法模拟GaAs衬底上InGaAs异质结表面形貌[J]. 材料导报2018,32(z1):547-552.

WANG G, WENG Y, WU P. Phase field simulation on the morphology of InGaAs heterostructure on GaAs substrate[J]. Materials Review, 2018, 32(z1): 547-552. (in Chinese)

[18] PROVATAS N, ELDER K. Phase-field methods in materials science and engineering[M]. Wiley-VCH Verlag GmbH & Co, 2010.

[19] BINER S B. Programming phase-field modeling[M]. Springer International Publishing, 2017.

[20] 徐慢, 夏冬林, 杨晟, 等．薄膜太阳能电池[J]．材料导报, 2006, 20(9): 109-111．

XU M, XIA D, YANG S, et al．Thin film solar cells[J]．Materials Review，2006, 20(9): 109-111．(in Chinese)

[21] 罗子江，周勋，杨再荣, 等．InGaAs/GaAs异质薄膜MBE生长研究[J]．功能材料, 2011, 42(5): 846-849.

LUO Z, ZHOU X, YANG Z, et al. The MBE growth research on InGaAs/GaAs heterofilms[J]. Journal of Functional Materials, 2011, 42(5): 846-849. (in Chinese)

[22] ARTHUR J R. Molecular beam epitaxy[J]. Surface Science, 2002, 500:189.

[23] PATRICK MCCRAY W. MBE deserves a place in the history books[J]. Nature Nanotechnology, 2007, 2: 259.

[24] BAUER E. Phanomenologische theorie der Kristallabscheidung an oberflachen. I[J]. Zeitschrift fur Kristallographie, 1958, 110: 372.

[25] STRANSKI I N, KRASTANOW L. Zur Theorie der orientierten Ausscheidung von Ionenkristallen aufeinander[J]. Abhandlungen der Mathematisch-Naturwissenschaftlichen Klasse IIb. Akademie der Wissenschaften Wien. 1938, 146: 797.

[26] MATTHEWS J W, BLAKESLEE A E. Defects in epitaxial multilayers-I: Misfit Dislocations[J]. Journal of Crystal Growth. 1974, 27: 118.

[27] PEOPLE R, BEAN J C. Calculation of crytical layer thickness versus lattice mismatch for GexSi1-x/Si Strained-Layer heterostructures[J]. Applied Physics Letters, 1985, 47: 322.

[28] SPECK S, POMPE W. Domain configurations due to multiple misfit relaxation mechanisms in epitaxial ferroelectric thin films. I. Theory[J]. Journal of Applied Physics,1994, 76: 466.

[29] MAREE P M J, BARBOUR J C, VAN DER VEEN J F, et al. Generation of misfit dislocations in semiconductors[J]. Journal of Applied Physics,1987, 62: 4413.

[30] SHENG G, HU J M, ZHANG J X, et al. Phase-field simulations of thickness-dependent domain stability in PbTiO3 thin films[J]. Acta Materialia, 2012, 60: 3296.

[31] SCHLOM D G, CHEN L Q, PAN X Q, et al. A thin film approach to engineering functionality into oxides[J]. Journal of the American Ceramic Society, 2008, 91(8): 2429.

[32] BEAN J C, FELDMAN L C, FIORY A T, et al. GexSi1-x/Si strained‐layer superlattice grown by molecular beam epitaxy[J]. Journal of Vacuum Science & Technology A: Vacuum, Surfaces, and Films, 1984, 2: 436.

[33] HOUGHTON D C, GIBBINGS C J, TUPPEN C G, et al. Equilibrium critical thickness for Si1-xGex strained layers on (100) Si[J]. Applied Physics Letters, 1990, 56, 460.

[34] ELMAN B, KOTELES E S, MALMAN P, ET AL. In situ measurements of critical layer thickness and optical studies of InGaAs quantum wells grown on GaAs substrates[J]. Applied Physics Letters, 1989, 55, 1659.

[35] AYERS J E. Heteroepitaxy of semiconductors: Theory, Growth, and Characterization[M]. CRC Press, 2017.

[36] ASARO R J, et al. Interface morphology development during stress corrosion cracking: Part I. Via surface diffusion[J]. Metallurgical Transactions, 1972, 3: 1789.

[37] GRINFELD M A. Instability of the equilibrium of a nonhydrostatically stressed body and a melt Soviet physics[J]. Doklady, 1986, 31: 831.

[38] VENABLES J A, SPILLER G D T, HANBUCKEN M. Nucleation and growth of thin films[J]. Reports on Progress in Physics, 1984, 47(4): 399.

[39] SROLOVITZ D J. On the stability of surfaces of stressed solids[J]. Acta Metallurgica, 1989, 37: 621.

[40] GAO H J. Some general properties of stress-driven surface evolution in a heteroepitaxial thin film structure[J]. Journal of the Mechanics and Physics of Solids, 1994, 42: 741.

[41] BALIBAR S, ALLES H, PARSHIN A Y. The surface of helium crystals[J]. Reviews of Modern Physics, 2005, 77, 317.

[42] JESSON D E, CHEN K M, PENNYCOOK S J, et al. Mechanism of strain induced roughening and dislocation multiplication in SixGe1-x thin films[J]. Journal of Electronic Materials, 1997, 26: 1039.

[43] BERREHAR J, CAROLI C, LAPERSONNE-MEYER C, et al. Surface patterns on single-crystal films under uniaxial stress: Experimental evidence for Grined instability[J]. Physical Review B, 1992, 46: 13487.

[44] KUZNETSOV P V, PANIN V E, PETRAKOVA I V. Grinfeld instability in the formation of a tweed structure at the Al crystal surface under cyclic tension[J]. Physical Mesomechanics, 2010, 13: 70.

[45] GURURAJAN M P, LAHIRI A. Elastic stress effects on microstructural instabilities[J]. Journal of the Indian Institute of Science, 2016, 96(3): 199.

[46] KHACHATURYAN A G. Theory of Structural Transfomations in Solid[M]. New York: Wiley, 1983.

[47] TAKAKI T, HIROUCHI T, TOMITA Y. Phase-field study of interface energy effect on quantum dot morphology[J]. Journal of Crystal Growth, 2008, 310: 2248.

[48] CHEN L Q, SHEN J. Applications of semi-implicit Fourier-spectral method to phase field equations[J]. Computer Physics Communications, 1998, 108: 147.

[49] CHUANG S L. Physics of optoelectronics devices[M]. New York: Wiley, 1995.

[50] MARIAGER S O, LAURIDSEN S L, DOHN A, et al. High-resolution three-dimensional reciprocal-space mapping of InAs nanowires[J]. Journal of Applied Crystallography, 2009, 42: 369.

[51] YU Y M, LIU B G. Phase-field model of island growth in epitaxy[J]. Physical Review E, 2004, 69: 021601.

[52] RATZ A, RIBALTA A, VOIGT A. Surface evolution of elastically stressed films under deposition by a diffuse interface model[J]. Journal of Computational Physics, 2006, 214: 187.

[53] ARJMAND M, DENG J, SWAMINATHAN N, et al. Effects of confinements on morphology of InxGa1-xAs thin film grown on sub-micron patterned GaAs substrate: Elastoplastic phase field model[J]. Journal of Applied Physics, 2014, 116: 114313.

[54] PODMANICZKY F, TóTH G I, TEGZE G, et al. Phase-field crystal modeling of heteroepitaxy and exotic modes of crystal nucleation[J]. Journal of Crystal Growth, 2017, 457: 24.

[55] WANG Y U, JIN Y M, KHACHATURYAN A G. Phase field microelasticity modeling of surface instability of heteroepitaxial thin films[J]. Acta Materialia, 2004, 52: 81.

[56] SEOL D J, HU S Y, LIU Z K, et al. Phase-field modeling of stress-induced surface instabilities in heteroepitaxial thin films[J]. Journal of Applied Physics, 2005, 98: 044910.

[57] YU P, HU S Y, CHEN L Q, et al. An iterative-perturbation scheme for treating inhomogeneous elasticity in phase-field models[J]. Journal of Computational Physics, 2005, 208: 34.

[58] NI Y, HE L H, SOH A K. Three-dimensional phase field simulation for surface roughening of heteroepitaxi-al films with elastic anisotropy[J]. Journal of Crystal Growth. 2005, 284: 281.

[59] LIANG X D, NI Y, HE L H. Shape-dependent composition profile in epitaxial alloy quantum dots: A phase-field simulation[J]. Computational Materials Science, 2010, 48: 871

[60] 大澤裕樹, 山中晃徳, Willy Kurnia, など. 超微細塑性加工と自己組織化による規則配列した金ナノドットアレイの創製[J]. 日本機械学会論文集(C編)，原著論文No.2010-JCR-0328

OSAWA H, YAMANAKA A, KURNIA W, et al. Fabrication of ordered gold nano dot array by nano plastic forming and self-assembly[J]. Transactions of the JSME (in Japanese), No.2010-JCR-0328

[61] ELDER K R, PROVATAS N, BERRY J, et al. Phase-field crystal modeling and classical density functional theory of freezing[J]. Physical Review B, 2007, 75: 064107.

[62] BOYNE A, DREGIA S A, WANG Y. Concurrent spinodal decomposition and surface roughening in thin solid films[J]. Applied Physics Letters, 2011, 99: 063111.

[63] BOYNE A, RAUSCHER M D, DREGIA S A, et al. Surface island formation with localized stresses[J]. Scripta Materialia, 2011, 64: 705.

[64] WANG L Y, LIU Z L, ZHUANG Z. Continuum modeling of surface roughening in heteroepitaxial structures based on phase field theory[J]. Computational Materials Science, 2017, 136: 109.

[65] WANG Y U, JIN Y M, KHACHATURYAN A G. Phase field microelasticity modeling of dislocation dynamics near free surface and in heteroepitaxial thin films[J]. Acta Materialia, 2003, 51: 4209.

[66] WANG Y U, JIN Y M, KHACHATURYAN A G. Mesoscale modelling of mobile crystal defects-dislocations, cracks and surface roughening: Phase field microelasticity approach[J]. Philosophical Magazine, 2005, 85: 261.

[67] WANG L Y, LIU Z L, ZHUANG Z. Developing micro-scale crystal plasticity model based on phase field theory for modeling dislocations in heteroepitaxial structures[J]. International Journal of Plasticity, 2016, 81: 267.

[68] HU S Y, CHEN L Q. Spinodal decomposition in a film with periodically distributed interfacial dislocations[J]. Acta Materialia, 2004, 52: 3069.

[69] YU P, HU S Y, CHEN L Q, et al. An iterative-perturbation scheme for treating inhomogeneous elasticity in phase-field models[J]. Journal of Computational Physics, 2005, 208: 34.

[70] HAATAJA M, MULLER J, RUTENBERG A D, et al. Dynamics of dislocations and surface instabilities in misfitting heteroepitaxial films[J]. Physical Review B, 2002: 035401.

[71] WU P P, GAO F L, ZHANG S G, et al. Surface morphology of GaAs/In0.3Ga0.7As in an elastic field of static point defects[J]. Chinese Physics Letters, 2014, 31(2): 026802.

[72] CHIRRANJEEVI B G, ABINANDANAN T A, GURURAJAN M P. A phase field study of morphological instabilities in multilayer thin films[J]. Acta Materialia, 2009, 57: 1060.

[73] ZAEEM M A, SINISA D J. Mesarovic morphological instabilities in thin films: Evolution maps[J]. Computational Materials Science, 2011, 50: 1030.

[74] WU P P, GAO F L, LI G Q. Effects of buffer layer thickness on the surface roughness of In0.3Ga0.7As thin films: A phase-field simulation[J]. Journal of Materials Research, 2013, 28(23).

[75] WU P P, WANG G. A phase-field model for multilayered heterostructure morphology[J]. Materials Science Forum. 2019, 944: 788-794.

[76] μ-pro package on Chen’s research group[EB/OL]. https://www.ems. psu.edu/～chen /index.html.[1] JOHN V. Introduction to surface and thin film processes[M]. Cambridge: Cambridge University Press, 2000.