类石墨烯氮化碳结构(C3N)热传导机理研究*

任国梁 申开波 刘永佳† 刘英光

1) (上海交通大学化学化工学院、材料科学与工程学院、分析测试中心,上海 200240)

2) (华北电力大学能源动力与机械工程学院,保定 071003)

类石墨烯氮化碳结构(C3N)作为一种全新的碳基二维半导体材料,由于其优异的机械和电子性能引起了研究者们的广泛关注,不同结构C3N 的热输运和声子输运机制还待进一步研究.本文构造了4 种不同结构的C3N,采用非平衡分子动力学与晶格动力学方法对不同结构的C3N 的热传导机理进行了研究.研究结果表明: 1)在4 种结构中M3 热导率最高,M1 次之,M4 热导率最低;2)不同结构的C3N 的热导率具有明显的尺寸效应和温度效应.当样本长度较短时,声子主要以弹道输运的方式进行传输;当样本长度增大,扩散输运占主导地位;随着温度的升高,Umklapp 散射在热输运中占据主导地位,使得热导率与温度具有1/T 的依赖性.3)与M3 相比,M1 和M4 结构中都存在更大的声子带隙,色散曲线进一步软化,低频和高频声子同时出现了局域化的特征,对热导率产生了显著的抑制作用.本文为更好地设计热管理材料提供了思路.

1 引 言

材料是人类文明的基石,人类的发展离不开对新材料的不断发现.自从2004 年曼彻斯特大学Novoselov 团队[1]成功分离出单层石墨烯,“碳”这种生命基础元素所构建出的二维材料成为了纳米材料研究的新领域,引发了前所未有的科研热潮.石墨烯具有蜂窝状的二维晶格结构,由于电子在石墨烯中可进行无质量狄拉克费米子行为[2],使得石墨烯具有极高的导电能力,大量的理论研究和实验数据都证明了石墨烯具有极高的载流子迁移能力(1010cm2/(V·s))[3−4],这让石墨烯在后硅CMOS 时代具备强大的竞争力.但是,石墨烯作为一种半金属材料并不存在带隙,这就意味着基于本征石墨烯的器件无开关性能[5].对于二进制逻辑运算来说,当微电子器件始终处于“开”或者“关”的状态时,便不具备利用价值[6].因此,通过对石墨烯的性能和空间结构进行调控,打开带隙,具有重大的科学意义.

为了解决这个问题,研究人员提出了诸如异质结构[7]、构建纳米网[8]以及缺陷工程[9]等研究方案.2007 年,以B,N 原子对石墨烯进行空穴和电子掺杂[10],得到了具有适合带隙的半导体性掺杂石墨烯,同时MoS、WS、硅烯、锗烯等非碳二维材料的相继出现也备受研究者的关注[11].但是以上技术手段的成本极其昂贵并且操作复杂,限制了此类技术的发展.因此碳基二维材料的改性研究重新得到重视和发展.

氮元素作为在元素周期表中与碳相邻的元素,具有与碳相似的原子半径,可以较为容易地与碳元素进行替代,形成多种新型二维层状结构.因此直接制备二维的类石墨烯结构便成为一个极具竞争力的方案.1996 年Teter 和Hemley[12]采用第一性原理对C3N4重新进行计算,提出C3N4具有5 种结构;2015 年,Mahmood 等[13]制备了一种多孔氮化石墨烯(C2N);2017 年,Yang 等[14]通过2,3−二氨基吩嗪在水热条件下的偶联反应,制备了无孔石墨烯状结构的新型二维半导体材料C3N,作为一种全新的碳基二维半导体材料,C3N 结构显示出比其他结构更高的热性能和机械性,其载流电子迁移率可达220 cm2/(V·s),并且基于单层C3N 薄膜的FET 器件开关比可以高达5.5 × 1010.Wei 等[15]通过控制堆叠顺序或施加电场来设计双层C3N 的带隙,实验结果表明双层C3N 可以提供可控的开/关比、高载流子迁移率和光电检测能力;同时,Mortazavi[16]通过分子动力学(molecular dynam−ics,MD)模拟出C3N 超晶格在300 K 的热导率为815 ± 20 W/(m·K),并且研究发现单层C3N 可以承受4000 K 的高温;Gao 等[17]使用基于第一性原理的玻尔兹曼输运方程(Boltzmann transport equation,BTE)计算发现C3N 的本征晶格热导率在室温下为380 W/(m·K);张浩课题组[18]预测单层C3N 室温下热导率为482 W/(m·K).但是在C3N中是否存在布里渊区折叠、声子局域化等一系列有趣的声子传输现象却有待进一步研究.

因此,本文构造了4 种不同结构的C3N,采用非平衡分子动力学法(non−equilibrium molecular dynamics,NEMD)研究C3N 材料的导热性能,探讨温度、结构等因素对其热导率的影响,并且利用晶格动力学对其声子输运过程进行研究,更加全面地分析C3N 结构的热输运机理.本文首先介绍了模型和模拟方法,包括计算热导率的分子动力学,然后给出了计算结果,计算了不同温度以及样本长度下的热导率,并详细讨论了其声子传输机理.

2 模型与方法

构造的4 种不同的C3N 结构如图1 所示.在C3N 的平面方向内,应用了周期性边界条件,同时为了避免周期性边界条件产生的层间相互作用,在Z方向上加入了20 Å的真空层.模型的最小重复单元包含4 × 4 个原始晶胞,由32 个原子组成.本文采用的纳米尺度传热试样的原理模型如图2所示,其中X方向为热流方向,为了防止与外界产生热量交换,在模型两端固定10 原子层厚度为固定边界,并且设定该区域粒子的速率为0.

图1 C3N 的模型结构示意图Fig.1.Schematic diagram of C3N model structure.

图2 NEMD 模拟计算热导率的原理图Fig.2.Schematic diagram of thermal conductivity calculated by NEMD simulation.

本次研究均是基于LAMMPS 软件进行的[19].在分子动力学的模拟中,势函数的选择决定了能否真实反应物质特性以及准确的描述原子(分子)间的相互作用力.因此采用Kinaci 等[20]改进以后的Tersoff 势函数,能更精准描述C3N 的物质特性以及共价键之间的相互作用力.描述C3N 性能的势参数(X 表示键相互作用的原子,此处X 可以是C 和N 中的任意一个)如表1 所示.

表1 C3N 原子相互作用的Tersoff 势函数参数[20]Table 1. Tersoff potential function parameters of C3N atom interactions.

在C3N 的NEMD 模拟中,将原子运动积分步长设置为0.5 fs,首先在0 K 的条件下对体系进行能量最小化,优化原子位置,使系统达到稳定状态,接着将系统置于等温等压系综(NPT)中弛豫3 ns,然后在正则系综(NVT)中弛豫3 ns,使系统结构松弛到平衡态.最后使用朗之万恒温器在体系中设置一个温度梯度,将体系放置于微正则系综(NVE)中进行6 ns.将体系沿着热流方向方向分为50 块,使用下式对模型内原子的平均温度进行统计:

其中,N为体系划分的每一部分原子数.随后将X方向上温度分布的线性部分进行拟合,得到温度梯度(∂T/∂x),如图3(a)所示.系统的热流密度J可由热源/冷源区域的能量交换率 dE/dt计算:

其中,A为模拟系统在垂直于热传输方向的横截面积,结果如图3(b)所示.

图3 C3N 中温度分布和能量分布 (a) 模型的温度分布;(b) 热源和冷源能量Fig.3.Temperature distribution and energy variation in C3N model: (a) Temperature distribution of the model;(b) heat source and heat sink energy.

最后,基于傅里叶导热定律计算体系热导率:

3 结果与分析

3.1 热导率

采用NEMD 方法计算不同结构的C3N 在300 K 下的热导率.计算结果如图4 所示.C3N 的长度被设定为25 nm,宽度保持为10 nm.根据Song 等[21]的研究,当样品宽度处于5—20 nm 范围内并采用周期性边界条件时,计算出的热导率的结果差异可以忽略不计,因此10 nm 宽度对于后续获得热输运的计算结果是合理的.从图4 可以看出,不同结构的C3N 之间热导率存在显著差异,完美结构(M3)的热导率最高,为368.77 W/(m·K),这与Kumar 等[22]模拟计算的结果基本一致.M1,M2,M4 的热导率分别为320.21 W/(m·K),200.70 W/(m·K),138.33 W/(m·K).从图1 可以明显看出,相比于M2,M1 和M3 两种结构拥有更大的颈部宽度,而颈宽的增大会显著提高材料热导率,因此M1 和M3 相比M2 拥有更高的热导率.在Wei 等[23]的研究中发现真实热流路径的伸长、热流面积的高估等将导致热导率降低,而M3 与M1 相比,M1 拥有更长的热流路径,故M1 结构的热导率的值也小于M3.在4 种结构中M4 的热导率最低,可推测这是由于在M4 中存在N—N 键,体系中键能差异较大,而在长波长的限制下,声子群速度决定了热导率,体系键能差异较大将增强声子散射,抑制声子热传导,从而降低热导率[21].

图4 不同结构的C3N 在300 K 下的热导率的热导率Fig.4.Thermal conductivity of C3N with different struc−tures at 300 K.

3.2 尺寸效应和温度效应

为了探究C3N 结构的热输运性能,计算了300 K 下不同样本长度下C3N 结构的热导率,结果如图5 所示.研究发现C3N 结构的热导率具有明显的尺寸效应[24],当模型样本长度范围在25—200 nm 时,热导率随着样本长度的增大急剧增大,这是因为在此过程中声子输运机制主要为弹道输运.随着样本长度的进一步增大,热导率增大的趋势降低,此时扩散输运占主导地位.当样本为800 nm 时,M1,M2,M3,M4 的热导率分别达到660.01 W/(m·K),467.95 W/(m·K),728.91 W/(m·K),371.94 W/(m·K).

图5 样本长度对热导率的影响Fig.5.Influence of sample length on thermal conductivity.

图6 显示了温度对C3N 热导率的影响,从图6可以看出,不同C3N 结构对温度的依赖性相似,表现出较强的非谐效应,这与An 等[25]的研究报道相似.根据声子的散射机制,温度的升高将会激发更多的高频声子,增大热导率,但同时Umklapp散射也会随之增强[26],从而降低声子群速度,导致热导率减小.因此对于C3N 结构,随着温度升高Umklapp 散射在热输运中占据主导地位,使得热导率与温度呈现1/T趋势.

图6 温度对热导率的影响Fig.6.Influence of temperature on thermal conductivity.

3.3 布里渊区折叠和声子局域化

为了进一步阐述热导率降低的机理,通过晶格动力学(general utility lattice program,GULP)[27]对不同C3N 结构中声子的振动本征模式进行分析.从图7 的色散曲线可以看出,不同结构的C3N 色散曲线存在较大差异.首先M1 和M3 相比,在红色椭圆突出的区域内M1 存在较大的带隙,而M3是连续的;M2 的色散曲线相对于M1 和M3 较为平坦,色散曲线的平坦意味着群速度的降低,对体系热传输起抑制作用;而M4 不仅存在较大的带隙(蓝色椭圆突出部分),声子带的折叠也致使色散曲线更为平坦,从而导致热导率更低[25].因此M3 结构更有利于热传输,M1 与M2 次之,M4 最差,这与之前的NEMD 计算得到的热导率结果相一致.

图7 C3N 的色散曲线Fig.7.Phonon dispersion curves for C3N.

为了定量分析热导率之间的差异,计算得到不同结构的声子群速度,如图8 所示.从图8 可以看出,不论高频光学支还是低频声学支,M4 的群速度都是最低的,在低频声学支,群速度更是从21.7 km/s 显著地降至17.1 km/s.同时在20—30 THz 之间,在图8 红色椭圆突显的区域内M1同M3 相比更大一部分的声子群速度集中在了5 km/s 以下,这极大程度上阻碍了声子的热传输.

图8 C3N 的群速度Fig.8.Phonon group velocities for C3N.

声子态密度(phonon density of states,PDOS)是一种用来表征其中声子活动的方法.通过对速度自关联函数(velocity autocorrelation function,VACF)进行傅里叶变换[28]计算得到的PDOS:

其中,N为原子总数;vi(t) 为第i个原子在t时刻的速度;vi(0) 为原子初速度.

图9 为计算得到的不同C3N 结构的声子态密度.可以看出C3N 结构在高频范围由平面内的纵向声学声子(LA 模声子)和横向声学声子(TA 模声子)主导,而平面外的声学声子(ZA 模声子)则主导低频.从图9 得知M1 面内的高频特征峰(G 峰)为0.038,M2 为0.036,M3 为0.044,M4 为0.029,相比与M3,其他3 个结构高频段的G 峰受到了不同程度地抑制,峰值急剧减小,其中M4 最明显,同时低频段的PDOS 峰也受到抑制,但其程度远小于高频波段.对于面内部分,M1 和M4 的G 峰由于平面内声子模式双重退化导致对称破坏分裂为两个峰[29];此外,M4 中的PDOS 中还出现了一些额外的小峰(图9 中蓝色箭头突出部分),这些额外的峰与声子局域化相对应[30].而造成声子局域化的原因是多重弹性散射的破坏性干扰,从而使得某些声子模式从传播转变为局域[31],这也是导致M4 热导率极低的主要原因之一.而在频率小于30 THz 的低频段,PDOS 贡献主要来自面外部分.通过对比图9 中面外部分的PDOS 可以看出,M3 在低频部分的所受到的抑制效果明显小于其他三者,这与图4 的结果相一致.

图9 C3N 的声子态密度Fig.9.Phonon density of states for C3N.

声子参与率是用以衡量参与某一特定模式的原子比例,在定量描述声子局域化效应得到很好的应用[32−33].为了进一步研究声子活动对不同结构C3N 热导率的影响,根据(6)式计算了不同结构的声子参与率(phonon participation ratio,PPR)[34]:

其 中,P DOSi(ω) 表示频率为ω的第i个原子的PDOS.

以PPR ＜ 0.4 时标准判断声子是否处于局域模,从图10 可看出,不管低频声子还是高频声子,M1 和M4 都表现出了较强局域化的特征,并且M4 中的局域化更为明显.这证明了M4 结构的会导致更大一部分声子产生局域化,而局域化的声子不能像离域声子那样有效的传递热量,因此导致M4 的热导率最低.而M2 在高频声子(尤其是50 THz)局域化程度较大,而这恰恰是PDOS 中特征G 峰对应的频率范围,高频活性较高的声子反而是处于局域的状态,这将显著降低材料的热导率.而M1 和M3 结构的参与率大部分都大于0.4,表明在M1 和M3 结构中声子模拟为传播的声子的比例较大,也更有利于声子热传输,这也与之前的计算结果相一致.

图10 C3N 的声子参与率Fig.10.Phonon participation ratios for C3N.

4 结 论

本文采用NEMD 方法研究了4 种不同结构的C3N 的热传输特性,系统地分析了不同结构的C3N 在相同周期长度下的热传导机理以及样本总长度和温度温度对热导率的影响,并详细对声子传输机理进行了讨论.结果表明:

1) 4 种结构中M3 热导率最高,M1 和M2 次之,M4 热导率最低.M1—M4 热导率的值分别为320.2 W/(m·K),200.7 W/(m·K),368.77 W/(m·K),138.33 W/(m·K).

2) C3N 具有明显的尺寸效应和温度效应.当样品尺寸较小时,随着模型长度的增大,热导率急剧增大,在此过程中声子主要以为弹道输运为主.当样品长度的进一步增大,由于扩散输运,热导率增大的趋势降低;Umklapp 散射在热输运中占据主导地位,导致C3N 结构热导率随温度的升高而降低.不同C3N 结构对温度的依赖性相似,表现出较强的非谐效应,并且随着温度升高Umklapp 散射在热输运中占据主导地位,使得热导率与温度呈1/T趋势.

3)相比于M3 结构,M1 存在较大的声子带隙,且由于平面内声子模式双重退化,导致对称破坏分裂为两个峰;并且不论低频声子还是高频声子都出现了较强局域化的特征,因此热导率小于M3.同M3 和M1 相比,M2 的特征G 峰被抑制,色散曲线较为平坦,群速度较小,抑制声子热传导;而M4 结构中也存在较大的带隙,色散曲线也更为平坦,不论高频光学支还是低频声学支,M4 群速度均最低,局域化程度更为明显,M4 结构热传输特性最差.