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

（湖南師范大學(xué)計算機部量子結(jié)構(gòu)與調(diào)控教育部重點實驗室，長沙410081）

1 引 言

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

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

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

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

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

3 溫度的量子化修正

在分子動力學(xué)模擬中，根據(jù)Boltzmann統(tǒng)計，分子動力學(xué)溫度可由以下計算公式得到：

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

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

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

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

定義速度自相關(guān)函數(shù)［8］：

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

結(jié)合（3）式及 （4）式，則有：

∞

根據(jù) （8）式，即可將分子動力學(xué)溫度TMD轉(zhuǎn)化為絕對溫度T，從而進行溫度修正.

4 結(jié)果分析與討論

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

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

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

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

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

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

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

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

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

［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］.大連：大連理工大學(xué)，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晶體薄膜構(gòu)件熱學(xué)性質(zhì)的晶格動力學(xué)研究［D］.長沙：國防科技大學(xué)，2008］