硅納米晶體薄膜熱膨脹性質的分子動力學研究

黃建平，楊 程，胡詩一，陽 波

（湖南師范大學計算機部量子結構與調控教育部重點實驗室，長沙410081）

1 引 言

硅納米晶體薄膜作為微納電子器件和微納機電系統中最常用的襯底材料，其熱膨脹性質的研究具有非常重要的意義：一方面硅納米晶體薄膜與其他薄膜熱膨脹系數的失配會產生熱應力，而較大的熱應力會引起系統的失效或損壞；另一方面，硅納米晶體薄膜的熱膨脹也能成為微熱執行器的動力來源.由于硅納米晶體薄膜的尺寸很小，對其熱膨脹性質進行實驗測量將非常困難.運用分子動力學 （MD），可以通過數值求解牛頓運動方程得到原子位置和速度等隨時間變化的信息，在此基礎上得到材料的熱學性質.分子動力學是溝通微觀與宏觀的橋梁，是研究熱學性質的一種非常有用的理論工具，目前已被成功地運用于納米晶體薄膜的熔化性質、比熱和熱傳導研究［1-3］，然而還未見過用于納米晶體薄膜的熱膨脹性質研究的報導.在硅納米晶體薄膜理想的 （100）表面，每個原子有兩個空鍵，空鍵的存在使表面相鄰的Si原子之間相互吸引，而有可能形成新的共價鍵，發生二聚，從而使表面重構［4-5］，并最終對其熱膨脹性質產生影響.本文將用分子動力學方法研究表面原子二聚對硅納米晶體薄膜熱膨脹性質的影響.

2 硅納米晶體薄膜的分子動力學模擬

我們選用了由美國Sandia國家實驗室開發的、開源免費的LAMMPS 軟件［6］，在八核計算機上對硅納米晶體薄膜進行分子動力學的并行模擬計算.

由于相對于其它勢函數，S-W 勢函數能夠較好地重現塊體硅單晶的熔點、熱傳導系數等熱學參數［7］，因此我們選擇了S-W 勢.硅納米晶體薄膜中原子的初始位置按其晶體結構給出，初始速度由設定溫度下的Maxwell-Boltzmann分布隨機選取；薄膜面向安排10×10個原胞，采用周期性邊界條件，薄膜法向安排10層原胞.系統為等溫等壓（NPT）系綜，采用恒零壓，對不同溫度下的硅納米晶體薄膜進行分子動力學模擬，時間步長為1.0×10-15s，先用1.0×106步進行弛豫以達到熱力學平衡，再用0.5×106步獲取所有硅原子的位置與速度隨時間的演化.通過在熱力學平衡狀態下對各個原子所有不同時刻的位置坐標值求平均，可以求得硅納米晶體薄膜中各個晶面的位置，并進而求得相應晶面間距，也能計算出晶格常數.硅納米晶體薄膜的熱膨脹系數可定義為：

其中L 表示硅納米晶體薄膜面向的晶格常數a或晶面間距d.

3 溫度的量子化修正

在分子動力學模擬中，根據Boltzmann統計，分子動力學溫度可由以下計算公式得到：

其中，＜…＞表示平均值，TMD表示分子動力學溫度，各個i表示不同的原子.由此得總的晶格振動內能與分子動力學溫度的關系為：

然而以上公式只有在量子效應可以忽略的情況下才成立，即要求溫度遠大于Debye溫度，在一般情況下，需要對分子動力學溫度進行量子化修正.

在固體量子理論中，晶格振動內能可表示為：

其中，ω 是聲子角頻率，ρ（ω）是聲子譜密度，n（ω，T）是聲子在絕對溫度T 下的布居數：

定義速度自相關函數［8］：

聲子譜密度可由速度自相關函數D（t）進行Fourier變換得到［8］：

結合（3）式及 （4）式，則有：

∞

根據 （8）式，即可將分子動力學溫度TMD轉化為絕對溫度T，從而進行溫度修正.

4 結果分析與討論

為了將硅晶體薄膜材料和塊體材料的熱膨脹性質進行對比，首先對其塊體材料進行MD 模擬.安排10×10×10個原胞，在三維方向上均采用周期性邊界條件，其它模擬參數與條件設置同納米晶體薄膜，模擬得到各不同溫度下塊體硅晶體的晶格常數，再對這些數據進行擬合，并對分子動力學溫度進行修正，得到塊體硅晶體的晶格常數隨溫度的變化關系曲線，有關結果如圖1所示.根據塊體硅晶體的晶格常數隨溫度的變化關系，得到如圖2所示的硅塊體材料熱膨脹系數隨溫度的變化關系，所得結果與Wei等人［9］的實驗測定值、Noya等人［10］基于Feynman路徑積分的MC 模擬結果以及黃建平等人［11］晶格動力學方法的計算結果吻合.這說明分子動力學模擬是能用于熱膨脹性質的研究并得到正確結果的.

圖1 硅塊體材料晶格常數隨溫度變化Fig.1 The lattice constant of Si bulk crystal vs temperature

圖2 硅塊體材料熱膨脹系數隨溫度變化Fig.2 The thermal expansion coefficient of Si bulk crystal vs temperature

同樣對硅納米晶體薄膜熱膨脹性質進行模擬，得到熱平衡狀態下其面向晶格常數和熱膨脹系數隨溫度的變化關系分別如圖3和圖4所示.由圖可知，與塊體材料相比，硅納米晶體薄膜熱的膨脹性質明顯不同，主要體現為：在約400K 以下的低溫段，隨著溫度的升高，硅納米晶體薄膜的面向的晶格常數反而減小，而相應的熱膨脹系數為負值.目前還未見文獻對該奇異現象有過報導.

一般地來講，熱膨脹是由原子之間的非和諧勢能引起的，在低溫時非和諧勢能較小，而在高溫時非和諧勢能較大，因此溫度的升高熱膨脹增加，即一般材料都呈現正的熱膨脹.納米硅晶體薄膜熱膨脹的奇異性質的存在肯定有基于其物質結構的內在原因.

圖3 硅納米晶體薄膜晶格常數隨溫度變化Fig.3 The lattice constant of Si nanocrystal film vs temperature

圖4 硅納米晶體薄膜熱膨脹系數隨溫度變化Fig.4 The thermal expansion coefficient of Si nanocrystal film vs temperature

納米硅晶體薄膜的每個內部硅原子的4 個SP3雜化軌道與周圍4個硅原子提供的4個SP3雜化軌道形成4個共價鍵，而納米硅晶體薄膜的每個表面硅原子只與次表面原子的形成2個共價鍵，另外在薄膜表面的外側還存在2個空鍵.在熱運動加劇的情況下，這些表面原子有可能在某個瞬間相互接近而使其空鍵對應的SP3雜化軌道有部分重疊，形成新的共價鍵，并使這些形成新的共價鍵的原子之間的距離縮短，從而使表面發生重構，這被稱為二聚現象［4-5］.塊體硅晶體的比表面積很小，其表面原子的二聚對熱膨脹性質的影響較小，可以忽略；而納米硅晶體薄膜的比表面積很大，其表面原子的二聚對熱膨脹性質的影響較大，必須加以考慮.

我們對處于不同溫度下的硅納米晶體薄膜表面層的二聚原子數目進行統計，發現在約400K 以下的低溫段，隨著溫度的升高薄膜表面層的二聚原子數目增加較多；而在約400K 以上的高溫段，由于表面原子中二聚的原子份額已經很大，再升高溫度雖然還會增加二聚原子數目，但增幅不會很大，并逐漸趨于飽和.于是在低溫下，隨著溫度的升高，由二聚現象而引起的晶格常數減小要比由原子之間的非和諧勢能引起晶格常數的增加要來得嚴重，從而導致硅納米晶體薄膜面向的晶格常數隨溫度增加而減小，故呈現負熱膨脹系數；而在高溫下，由原子之間的非和諧勢能引起的晶格常數的增加要比由二聚現象而引起的晶格常數減小要來得多，從而導致硅納米晶體薄膜面向的晶格常數隨溫度增加而增大，故整體呈現正熱膨脹系數.

［1］ Yang X，Wu D.The melting behaviors of the Nb（1 1 0）nanofilm：a molecular dynamics study［J］.Appl.Surf.Sci.，2010，256（10）：3197.

［2］ Adiga S，Adiga V P，Carpick R W，et al.Vibrational properties and specific heat of ultrananocrystalline diamond：molecular dynamics simulations［J］.J.Phys.Chem.C，2011，115（44）：21691.

［3］ Ju S，Liang X.An atomic level study on the out-ofplane thermal conductivity of polycrystalline argon nanofilm［J］.Chin.Sci.Bull.，2012，57：294.

［4］ Battaglia C，Gaál-Nagy K，Monney C，et al.Elementary structural building blocks encountered in silicon surface reconstructions［J］.J.Phys.：Condens.Matter，2009，21（1）：013001.

［5］ Fauster T，Tanaka S，Tanimura K.Unoccupied dimer-bond state at Si（001）surfaces［J］.Phys.Rev.B，2011，84：235444.

［6］ Plimpton S.Fast parallel algorithms for short-range molecular dynamics［J］.J.Comp.Phys.，1995，117：1.

［7］ Stillinger F H，Weber T A.Computer simulation of local order in condensed phases of silicon［J］.Phys.Rev.B，1985，31：5262.

［8］ Yu J.Fabrication of Microcalorimeter and Investigation on the Specific Heat of Micro／nanometer Thin Films［D］.Dalian：Dalian University of Technology，2006（in Chinese）［余雋.微量熱計的研制及微納米薄膜比熱測試與分析［D］.大連：大連理工大學，2006］

［9］ Wei S，Li C，Chou M Y.Ab initio calculation of thermodynamic properties of silicon［J］.Phys.Rev.B，1994，50：14587.

［10］ Noya J C，Herrero C P，Ramirze R.Thermodynamic properties of c-Si derived by quantum path-integral Monte Carlo simulations［J］.Phys.Rev.B，1996，53：9869.

［11］ Huang J P.Lattice dynamics study on thermal properties of crystal film in MEMS ［D］.Changsha：National University of Defense Technology，2008（in Chinese）［黃建平.MEMS晶體薄膜構件熱學性質的晶格動力學研究［D］.長沙：國防科技大學，2008］