掺杂硅纳米薄膜法向热导率的分子动力学模拟

陈 慧 魏志勇 陈伟宇 刘晨晗 毕可东 陈云飞

(1东南大学机械工程学院, 南京 211189)(2江苏海事职业技术学院船舶与海洋工程学院, 南京 211170)

掺杂硅纳米薄膜法向热导率的分子动力学模拟

陈 慧1,2魏志勇1陈伟宇1刘晨晗1毕可东1陈云飞1

(1东南大学机械工程学院, 南京 211189)(2江苏海事职业技术学院船舶与海洋工程学院, 南京 211170)

采用非平衡态分子动力学方法计算了平均温度为300 K时膜厚2.2～104.4 nm的硅纳米薄膜以及掺锗的硅纳米薄膜的法向热导率.采用随机掺杂方式在硅纳米薄膜中掺入锗原子,掺杂浓度分别为5%和50%.模拟结果表明,相同膜厚的掺锗硅薄膜的法向热导率比纯硅薄膜的法向热导率小很多,掺杂前后的硅薄膜法向热导率均随着膜厚的增大而增大.通过计算体态硅热导率关于声子平均自由程的累积函数,发现对体态硅热导率主要贡献是声子平均自由程为2～2 000 nm的声子,而掺锗的体态硅中对热导率贡献约占80%的声子平均自由程为0.1～30 nm,远小于纯硅中对热导率主要贡献的声子平均自由程.

硅纳米薄膜;热导率;掺杂;分子动力学;声子平均自由程;热导率累积函数

微电子器件的日益小型化使得其特征尺寸已缩小至纳米量级.一方面,对高密度热量的有效转移已成为影响微电子器件工作性能和可靠性的关键因素,因而期望构成微电子器件的纳米材料具有较高的热导率;另一方面,提升热电器件制冷效率的一个重要途径是在不影响热电材料电导率和Seebeck系数的情况下降低热导率,将显著提高热电材料的品质因数[1-2],大幅提升热电器件的性能.目前热电设备中使用较多的硅材料,其热电性能较差.因此,有效降低硅材料的热导率以提高其品质因数一直是该领域的研究热点之一.研究发现,当纳米器件或结构的特征尺寸降至与材料中传递热量的声子平均自由程(MFP)相当时,纳米结构的热传导将出现显著的尺寸效应[3].此时,材料的热导率主要取决于结构尺寸.由于电子的MFP比声子的MFP小得多,当材料尺度接近声子的MFP时,声子热导率会降低,而电子的电导率则影响较小.因此,利用纳米结构对声子输运的影响,降低材料的热导率,可有效提高热电材料的能量转换效率.纳米材料导热出现尺寸效应的主要原因在于,当材料的特征尺寸远小于声子平均自由程时,基本不会发生声子-声子间散射,导热处于弹道输运阶段,热导率随着尺寸的增加而增大.例如,Feng等[4]采用非平衡态分子动力学(NEMD)方法模拟了硅薄膜法向热导率,发现在膜厚为2～10 nm范围内,硅薄膜热导率显著低于体态实验值,且随着薄膜厚度的增大而增大.Gomes等[5]采用平衡态分子动力学(EMD)方法模拟了薄膜厚度为2～217 nm的硅薄膜面向和法向热导率.苏高辉等[6]采用蒙特卡洛方法研究了硅薄膜法向热导率出现明显尺度效应的临界尺寸,当薄膜厚度等于1 μm 时,开始呈现明显的尺度效应,此时扩散输运和弹道输运共存.Sun等[7]采用分子动力学方法研究了硅中声子输运的尺寸效应,认为硅薄膜热导率下降的主要原因是由于声子的限制和界面热阻的影响.Yang等[8]采用NEMD方法研究了硅纳米结构热导率的尺寸效应,发现硅纳米结构的长度对其热导率具有显著影响,并将其归因于硅纳米结构中声子有效平均自由程的显著减小以及纳米结构尺寸的限制导致对热传导作贡献的声子数量的减少.此外,硅的晶格缺陷以及掺杂等因素也会增加对声子的散射几率而使声子平均自由程减小,从而降低材料的热导率.Asheghi等[9]采用实验方法测量了温度为15～300 K时分别掺有硼和磷的单晶硅薄膜的热导率,结果表明杂质的存在显著降低了材料的热导率.Bi等[10]采用分子动力学模拟方法研究了在氩晶体中掺杂氪后对其热导率的影响,结果表明氩晶体的热导率降低了近1.9倍.Garg等[11]的研究表明硅和锗之间的质量差引起的质量失配对于抑制硅锗合金中的热传导起到了重要作用.目前关于掺锗的硅纳米薄膜的导热性能研究较少,且对热导率做贡献的声子平均自由程的大小尚未获得一致的结论.

本文采用NEMD方法研究了掺锗对硅纳米薄膜热导率的影响,讨论了热导率随膜厚的变化关系,计算了热导率关于声子平均自由程的累积函数[12],得出了对热导率做贡献的声子平均自由程的范围.

1 分子动力学模型

在硅纳米薄膜中以随机方式掺入锗(Ge)原子,原子的初始位置按照金刚石晶格结构排列在格点上,采用SW作用势[13]描述硅原子之间、锗原子之间以及硅原子与锗原子之间的相互作用力.图1为随机掺锗的硅纳米薄膜法向热导率的NEMD模型示意图,模拟区域大小用Nxa×Nya×Nza表示,其中Nx,Ny和Nz分别为x,y和z方向上的晶胞(unit cell, UC)数,a为硅的晶格常数.x方向为热流方向,分别在y和z方向上施加周期性边界条件,以实现对无限大平面纳米薄膜的模拟.模拟发现,垂直于热流方向的横截面积大小对模拟结果影响很小,本文中模拟区域横截面积均为10 UC×10 UC.在模拟区域的两端分别设置厚度为2 UC的固定绝热墙,热、冷域的厚度也均为2 UC,模拟样品的膜厚均不计算固定墙和热、冷域的厚度.

图1 掺杂硅纳米薄膜法向热导率的NEMD模型

模拟系统内原子速度的初始化遵循Maxwell分布,采用Velocity-Verlet算法对运动方程进行时间积分,积分步长Δt=0.5 fs,整个MD过程共模拟1.05×107步(5.25 ns)以上.最初的2.0×106步使模拟系统在NVT系综下演进,通过调整原子的速度使模拟系统的温度为恒定值.随后,系统在NVE系综下运行5.0×105步,保持系统能量守恒.此后至少再模拟8.0×106步,通过加减能量法在模拟区域内形成温度梯度,即每次给热域内增加能量Δε,同时从冷域内减少能量Δε.系统达到稳定状态时再运行6.0×106步以统计薄膜法向(x方向)的温度梯度.将整个模拟区域沿x方向划分为2Nx层,根据能量均分定理可计算每一层的区域温度Tj为

(1)

式中,m为原子质量;Nj为第j层所含的原子数;vi为原子i的速度;kB为Boltzmann常数.每隔1 000步统计一次各模拟层的温度,将最后运行的6.0×106步内统计的温度取平均值作为各层在热平衡后的平均温度,温度分布曲线如图2所示.根据分布曲线可计算出温度梯度▽T.模拟区域的法向热流密度Jx为

(2)

式中,A为薄膜的横截面积;Δt为积分步长.依据Fourier定律,掺锗硅纳米薄膜的法向热导率k为

(3)

图2 掺锗硅纳米薄膜的温度分布曲线

2 模拟结果与分析

本文分别模拟了300 K时锗原子的掺杂浓度为5%(Si0.95Ge0.05)、50%(Si0.5Ge0.5)以及不掺杂时的硅纳米薄膜在不同膜厚(2.2～104.4 nm)时的法向热导率,模拟结果如图3所示.在相同膜厚时,掺锗硅薄膜的法向热导率均显著低于纯硅薄膜的法向热导率.Yang等[14]采用NEMD方法对同位素掺杂的硅纳米线的热导率研究表明,随掺杂浓度的增加,热导率呈指数变化,极小的掺杂浓度即会引起热导率显著降低,而当热导率下降到一个最低值后,会随着掺杂浓度的增大而增大.

图3 硅纳米薄膜法向热导率随膜厚的变化(T=300 K)

硅薄膜中对热传导的贡献主要来自于声子的输运,忽略自由电子的作用,薄膜内除了声子与声子间散射外,薄膜边界以及晶格缺陷和杂质等对声子的散射占主导地位.在掺锗硅薄膜中,由于锗原子和硅原子质量以及直径的差异引发了声子-杂质间散射,而且,随机分布的大质量原子增加了硅薄膜中的晶格缺陷,使声子群速率降低,声子自由度减小,因此,掺锗大大降低了硅薄膜的热导率.

从图3可看出,当温度为300 K时,2种掺锗浓度以及不掺锗的硅纳米薄膜的法向热导率均随着膜厚的增加而增大,呈现出明显的尺寸效应.

如图4所示,采用基于Matthiessen的倒数拟合方法[15]外推,有1/k=A/t+B,其中,k为热导率，t为硅薄膜的法向厚度,当1/t=0时,1/B即为体态硅的热导率kbulk,A/B为声子平均自由程Λbulk,具体数据见表1.当温度为300 K时,体态硅的热导率为153.14 W/(m·K),声子平均自由程为110.7 nm.掺锗浓度为5%和50%的体态硅热导率分别为5.84和1.62 W/(m·K).

关于硅的声子平均自由程有多种不同的研究结果.根据声子气动理论计算出室温下硅晶体的声子平均自由程为41 nm[16];Alvarez等[17]根据传统的Debye模型,计算出室温时硅的声子有效平均自由程大约为40～43 nm;而Dames等[18]研究认为,硅的声子平均自由程为210 nm左右.为了明确不同平均自由程的声子对热导率的贡献,最近研究者[12, 19-21]提出了热导率关于声子平均自由程的累积函数计算方法，并进行了深入研究.

图4 硅纳米薄膜法向热导率倒数拟合(T=300 K)

物理量不掺锗掺杂锗浓度/%550kbulk/(W·(m·K)-1)153.145.841.62Λbulk/nm110.708.8010.50

(4)

式中,Λ为体态硅中三声子散射平均自由程.根据文献[19]可计算出

(5)

从而得到需求解的Fredholm积分方程为

(6)

式中,kfilm(t)为膜厚为t的硅薄膜法向热导率,由NEMD模拟获得;α(Λ)为热导率关于声子平均自由程的累积函数,表示声子平均自由程小于Λ的声子对热导率贡献的百分数.采用类似文献[20]的方法,先给定一个初值,然后采用一种优化方法求解方程(6)[12],得到

(7)

式中,Λ0可用表1中的Λbulk代替.图5给出了300 K时不掺锗、掺锗浓度为5%和50%的体态硅的热导率累积函数与声子平均自由程的关系.

由图5可看出,对于不掺锗的体态硅,热导率主要取决于声子平均自由程为2～2 000 nm的声子,其对热导率的贡献约占80%,声子平均自由程大于40 nm的声子对热导率的贡献超过50%,这与Regner等[21]采用测量方法获得晶体硅的热导率累积函数的研究结果较为一致.由图5可知,声子平均自由程大于100 nm的声子对热导率的贡献仍约占40%.因此,当硅薄膜的模拟厚度小于100 nm时,大多数声子处于弹道热输运机制,薄膜边界对声子的散射起主导作用,声子的有效平均自由程与膜厚相关,导致硅薄膜法向热导率的NEMD模拟值呈现明显的尺寸效应.

图5 热导率对声子平均自由程的累积函数

对于掺锗的体态硅,掺锗浓度为5%和50%的热导率累积函数随声子平均自由程的变化曲线非常接近.图5表明,掺锗的体态硅中对热导率的贡献80%均来自于声子平均自由程为0.1～30 nm的声子.对掺锗硅热导率作主要贡献的声子平均自由程远小于不掺锗硅的声子平均自由程,这是由于声子-杂质间的散射引起的.

3 结论

1) 采用NEMD方法研究了300 K时硅纳米薄膜以及掺锗的硅纳米薄膜的法向热导率随膜厚的变化关系，相同膜厚的掺锗硅薄膜的法向热导率比纯硅薄膜的小很多.

2) 不论掺锗与否，硅薄膜的法向热导率均随着膜厚的增大而增大，呈现出明显的尺寸效应.

3) 采用倒数拟合方法得到体态硅的热导率为153.14 W/(m·K)，掺锗浓度为5%和50%的体态硅热导率分别为5.84和1.62 W/(m·K).

4) 通过计算热导率关于声子平均自由程的累积函数，分析了对热导率做主要贡献的声子平均自由程.当温度为300 K时，对体态硅热导率作主要贡献来自于声子平均自由程为2～2 000 nm的声子，声子平均自由程大于100 nm的声子对热导率的贡献约占40%.这表明，传统观点严重低估了对体态硅热导率贡献的声子平均自由程的大小.掺锗浓度为5%和50%时，对热导率贡献约占80%的声子平均自由程均为0.1～30 nm，远小于纯硅中对热导率作主要贡献的声子平均自由程.

References)

[1]Liu W S, Yan X, Chen G, et al. Recent advances in thermoelectric nanocomposites[J].NanoEnergy, 2012, 1(1): 42-56. DOI:10.1016/j.nanoen.2011.10.001.

[2]Hochbaum A I, Chen R K, Delgado R D, et al. Enhanced thermoelectric performance of rough silicon nanowires[J].Nature, 2008, 451(7175): 163-167. DOI:10.1038/nature06381.

[3]Chen G.Nanoscaleenergytransportandconversion[M]. Oxford,UK: Oxford University Press, 2005: 283-291.

[4]Feng X L, Li Z X,Guo Z Y. Molecular dynamics simulation of thermal conductivity of nanoscale thin silicon films[J].MicroscaleThermophysicalEngineering, 2003, 7(2): 153-161. DOI:10.1080/10893950390203332.

[5]Gomes C J, Madrid M, Goicochea J V, et al. In-plane and out-of-plane thermal conductivity of silicon thin films predicted by molecular dynamics[J].JournalofHeatTransfer, 2006, 128(11): 1114-1121. DOI:10.1115/1.2352781.

[6]苏高辉,杨自春,孙丰瑞. 硅薄膜导热系数微尺度效应的临界尺寸[J]. 纳米技术与精密工程, 2014,12 (3): 176-181. DOI:10.13494/j.npe.20130086. Su Gaohui, Yang Zichun, Sun Fengrui. Critical size for microscale effect of silicon film thermal conductivity[J].NanotechnologyandPrecisionEngineering, 2014, 12(3): 176-181. DOI:10.13494/j.npe.20130086. (in Chinese)

[7]Sun L, Murthy J Y. Domain size effects in molecular dynamics simulation of phonon transport in silicon[J].AppliedPhysicsLetters, 2006, 89(17): 171919. DOI:10.1063/1.2364062.

[8]Yang Y W, Liu X J,Yang J P. Nonequilibrium molecular dynamics simulation for size effects on thermal conductivity of Si nanostructures[J].MolecularSimulation, 2008, 34(1): 51-55. DOI:10.1080/08927020701730419.

[9]Asheghi M, Kurabayashi K, Kasnavi R, et al. Thermal conduction in doped single-crystal silicon films[J].JournalofAppliedPhysics, 2002, 91(8): 5079-5088. DOI:10.1063/1.1458057.

[10]Bi K D, Zhao Y Y, Chen Y F, et al. The effects of different doping patterns on the lattice thermal conductivity of solid Ar[J].JournalofPhysicsandChemistryofSolids, 2012, 73(2): 204-208. DOI:10.1016/j.jpcs.2011.11.001.

[11]Garg J, Bonini N, Kozinsky B, et al. Role of disorder and anharmonicity in the thermal conductivity of silicon-germanium alloys: A first-principles study[J].PhysicalReviewLetters, 2011, 106(4): 045901. DOI:10.1103/PhysRevLett.106.045901.

[12]Wei Z, Yang J, Chen W, et al. Phonon mean free path of graphite along the c-axis[J].AppliedPhysicsLetters, 2014, 104(8): 081903. DOI:10.1063/1.4866416.

[13]Stillinger F H,Weber T A. Computer simulation of local order in condensed phases of silicon[J].PhysicalReviewB, 1985, 31(8): 5262-5271. DOI:10.1103/physrevb.31.5262.

[14]Yang N, Zhang G, Li B W. Ultralow thermal conductivity of isotope-doped silicon nanowires[J].NanoLetters, 2008, 8(1): 276-280. DOI:10.1021/nl0725998.

[15]Schelling P K, Phillpot S R,Keblinski P. Comparison of atomic-level simulation methods for computing thermal conductivity[J].PhysicalReviewB, 2002, 65(14): 144306. DOI:10.1103/physrevb.65.144306.

[16]Henry A S,Chen G. Spectral phonon transport properties of silicon based on molecular dynamics simulations and lattice dynamics[J].JournalofComputationalandTheoreticalNanoscience, 2008, 5(2): 141-152. DOI:10.1166/jctn.2008.2454.

[17]Alvarez F X,Jou D. Memory and nonlocal effects in heat transport: From diffusive to ballistic regimes[J].AppliedPhysicsLetters, 2007, 90(8): 083109. DOI:10.1063/1.2645110.

[18]Dames C, Chen G. Theoretical phonon thermal conductivity of Si/Ge superlattice nanowires[J].JournalofAppliedPhysics, 2004, 95(2): 682-693.

[19]Yang F,Dames C. Mean free path spectra as a tool to understand thermal conductivity in bulk and nanostructures[J].PhysicalReviewB, 2013, 87(3): 035437. DOI:10.1103/physrevb.87.035437.

[20]Minnich A J. Determining phonon mean free paths from observations of quasiballistic thermal transport[J].PhysicalReviewLetters, 2012, 109(20): 205901. DOI:10.1103/PhysRevLett.109.205901.

[21]Regner K T, Sellan D P, Su Z, et al. Broadband phonon mean free path contributions to thermal conductivity measured using frequency domain thermoreflectance[J].NatureCommunications, 2013, 4: 1640-01-1640-07. DOI:10.1038/ncomms2630.

[22]Majumdar A. Microscale heat-conduction in dielectric thin-films[J].JournalofHeatTransfer, 1993, 115(1): 7-16. DOI:10.1115/1.2910673.

Molecular dynamics simulation of out-of-plane thermal conductivity of doped silicon nanofilms

Chen Hui1, 2Wei Zhiyong1Chen Weiyu1Liu Chenhan1Bi Kedong1Chen Yunfei1

(1School of Mechanical Engineering, Southeast University, Nanjing 211189, China)(2Department of Ship and Ocean Engineering, Jiangsu Maritime Institute, Nanjing 211170, China)

The out-of-plane thermal conductivities of both pure silicon nanofilm and that doped with germanium (Ge) with film thicknesses ranging from 2.2 to 104.4 nm were calculated by non-equilibrium molecular dynamics (NEMD) simulation at an average temperature of 300 K. Ge atoms were doped in a random pattern with a doping density of 5% and 50%, respectively. The results show that there is a large reduction in thermal conductivities of the silicon nanofilm after doping with the same thickness. The out-of-plane thermal conductivity of the Ge-doping silicon nanofilm as well as that of a pure one increases with the increase of the film thickness. By calculating the thermal conductivity accumulation function in terms of phonon mean free path (MFP), it indicates that at 300 K the dominant contribution to the bulk thermal conductivity of pure silicon comes from phonons with MFPs between 2 and 2 000 nm, while the MFPs of phonons that contribute about 80% to the bulk thermal conductivity of Ge-doping silicon are between 0.1 and 30 nm, thus they are much less than those of the pure silicon.

silicon nanofilm; thermal conductivity; doping; molecular dynamics; phonon mean free path; thermal conductivity accumulation function

10.3969/j.issn.1001-0505.2017.03.013

2016-09-16. 作者简介: 陈慧(1982—)，女，博士生；陈云飞(联系人)，男，博士，教授，博士生导师，yunfeichen@seu.edu.cn.

国家自然科学基金资助项目(51205061, 51375092).

陈慧，魏志勇，陈伟宇，等.掺杂硅纳米薄膜法向热导率的分子动力学模拟[J].东南大学学报(自然科学版),2017,47(3):490-494.

10.3969/j.issn.1001-0505.2017.03.013.

TK124

A

1001-0505(2017)03-0490-05