基于准简谐格林-久保理论结合流体动力学外推方法的非晶二氧化铪导热机制研究

上一篇

下一篇

朱雄飞, 孙健时, 熊玉成, 李寿航, 刘向军. 基于准简谐格林-久保理论结合流体动力学外推方法的非晶二氧化铪导热机制研究[J]. 物理学报, 2025, 74(11): 116302-1. doi: 10.7498/aps.74.20250350
引用本文: 朱雄飞, 孙健时, 熊玉成, 李寿航, 刘向军. 基于准简谐格林-久保理论结合流体动力学外推方法的非晶二氧化铪导热机制研究[J]. 物理学报, 2025, 74(11): 116302-1. doi: 10.7498/aps.74.20250350
Xiongfei ZHU, Jianshi SUN, Yucheng XIONG, Shouhang LI, Xiangjun LIU. Research on thermal transport mechanism of amorphous hafnia based on quasi-harmonic Green-Kubo theory combined with hydrodynamic extrapolation method[J]. Acta Physica Sinica, 2025, 74(11): 116302-1. doi: 10.7498/aps.74.20250350
Citation: Xiongfei ZHU, Jianshi SUN, Yucheng XIONG, Shouhang LI, Xiangjun LIU. Research on thermal transport mechanism of amorphous hafnia based on quasi-harmonic Green-Kubo theory combined with hydrodynamic extrapolation method[J]. Acta Physica Sinica, 2025, 74(11): 116302-1. doi: 10.7498/aps.74.20250350

基于准简谐格林-久保理论结合流体动力学外推方法的非晶二氧化铪导热机制研究

    通讯作者: E-mail: shouhang.li@universite-paris-saclay.fr.;  E-mail: xjliu@dhu.edu.cn.
  • 中图分类号: 63.50.Lm, 66.70.-f, 02.70.Ns

Research on thermal transport mechanism of amorphous hafnia based on quasi-harmonic Green-Kubo theory combined with hydrodynamic extrapolation method

    Corresponding authors: E-mail: shouhang.li@universite-paris-saclay.fr.;  E-mail: xjliu@dhu.edu.cn.
  • MSC: 63.50.Lm, 66.70.-f, 02.70.Ns

  • 摘要: 非晶态材料二氧化铪在微电子器件中具有广泛应用, 理解其微观导热机制对于提升电子器件的性能和可靠性至关重要. 以往的研究大多基于分子动力学和单一准简谐格林-久保方法, 难以准确考虑低频振动模式的导热贡献. 本文基于准简谐格林-久保理论, 结合流体动力学外推法, 对不同有序度的非晶二氧化铪结构的热输运机制进行全面研究. 该方法可有效克服单一准简谐格林-久保方法中的有限尺寸问题. 理论预测表明, 非晶二氧化铪的热导率与微观结构有序度呈现弱相关性. 基于模态分析表明, 中低频振动模式对热导率具有显著贡献, 是单一准简谐格林-久保方法低估非晶二氧化铪热导率的主要原因. 同时, 本文基于非谐动态结构因子分离了传播子和扩散子对非晶二氧化铪导热的贡献, 计算表明扩散子在所有非晶二氧化铪结构导热中均占据主导作用. 然而, 传播子的导热贡献仍不可忽略, 其占比可高达20%以上, 且随着有序度的增大而增大.
  • 加载中
  • 图 1  生成非晶结构过程中势能和均方位移累计值随模拟时间的变化, 图中灰色虚线之前为熔融过程, 灰色虚线与黑色实线之间为在温度4000 K下的弛豫过程, 红色点划线之后为退火过程

    Figure 1.  Variation of potential energy and accumulated mean square displacement with simulation time during the generation of amorphous structures, the melting process is before the gray dashed line, the relaxation process at 4000 K is between the gray dashed line and the black solid line, and the annealing process is after the red dotted line.

    图 2  非晶HfO2在不同淬火速度下的(a)径向分布函数和(b)环数分布

    Figure 2.  (a) Radial distribution function and (b) ring distribution of amorphous HfO2 at different quenching rates.

    图 3  有限尺寸(2592个原子)和流体动力学外推法的非晶HfO2随温度变化的热导率, 图中图例的数值表示淬火速度, 误差棒为热导率标准差

    Figure 3.  Thermal conductivity of amorphous HfO2 with temperature for finite size (2592 atoms) and the hydrodynamic extrapolation. The values of the legend denote the quenching rate, error bar is standard deviations of thermal conductivity.

    图 4  非晶HfO2热导率随温度变化关系, 图中散点是Panzer等[45]、Lee等[46]、Scott等[10]和Chaubey等[47]报道的非晶HfO2实验数据

    Figure 4.  Temperature-dependent thermal conductivity of amorphous HfO2, scattered points in the figure are experimental data for amorphous HfO2 reported by Panzer et al.[45], Lee et al.[46], Scott et al.[10] and Chaubey et al.[47].

    图 5  300 K温度下不同淬火速度非晶HfO2 (a) 热容量; (b) 弛豫时间; (c) 模态扩散率; (d) 参与比

    Figure 5.  (a) Heat capacity; (b) lifetime; (c) modal diffusivity; (d) participation ratio of amorphous HfO2 with different quenching rates at 300 K.

    图 6  (a)—(e) 淬火速度分别为5×1012, 1×1012, 5×1011, 1×1011, 5×1010 K/s的非晶HfO2的纵向非谐动态结构因子; (f)—(j) 为对应淬火速度下的横向非谐动态结构因子; 图中颜色条是根据(8)式计算的非谐动态结构因子强度, 蓝色直线是非谐动态结构因子在低频部分的线性拟合, 红色虚线是选取的截止频率

    Figure 6.  (a)–(e) Longitudinal anharmonic dynamic structure factors of amorphous HfO2 with different quenching rates of 5×1012, 1×1012, 5×1011, 1×1011, and 5×1010 K/s; (f)–(j) transverse anharmonic dynamic structure factors at corresponding quenching rates. The colorbar indicates the strength of the anharmonic dynamic structure factor calculated by Eq. (8), the blue straight line is a linear fit of the anharmonic dynamic structure factor in the low frequency part, and the red dashed line is the selected cut-off frequency.

    图 7  300 K温度下不同淬火速度非晶HfO2的传播子和扩散子对热导率的贡献, 图中的蓝色文字代表传播子对总热导率的贡献百分比

    Figure 7.  Contributions of propagon and diffuson to thermal conductivity of amorphous HfO2 with different quenching rates at 300 K, the blue text in the figure represents the percentage contribution of propagon to the total thermal conductivity.

  • [1] Wilk G D, Wallace R M, Anthony J 2001 J. Appl. Phys. 89 5243 doi: 10.1063/1.1361065
    [2] Chen J H, Liang M F, Song Y, Yuan J J, Zhang M Y, Luo Y M, Wang N N 2024 Chin. Phys. B 33 047503 doi: 10.1088/1674-1056/ad1a88
    [3] Wang Y, Zahid F, Wang J, Guo H 2012 Phys. Rev. B 85 224110 doi: 10.1103/PhysRevB.85.224110
    [4] 董典萌, 汪成, 张清怡, 张涛, 杨永涛, 夏翰驰, 王月晖, 吴真平 2023 物理学报 72 097302 doi: 10.7498/aps.72.20222222 Dong D M, Wang C, Zhang Q Y, Zhang T, Yang Y T, Xia H C, Wang Y H, Wu Z P 2023 Acta. Phys. Sin. 72 097302 doi: 10.7498/aps.72.20222222
    [5] Gao R, Liu C, Shi B, Li Y, Luo B, Chen R, Ouyang W, Gao H, Hu S, Wang Y 2024 Chin. Phys. Lett. 41 087701 doi: 10.1088/0256-307X/41/8/087701
    [6] Chen S, Chen M, Liu Y, Cao D, Chen G 2024 Chin. Phys. B 33 098701 doi: 10.1088/1674-1056/ad4ff4
    [7] McKenna K, Shluger A, Iglesias V, Porti M, Nafría M, Lanza M, Bersuker G 2011 Microelectron. Eng. 88 1272 doi: 10.1016/j.mee.2011.03.024
    [8] 袁国亮, 王琛皓, 唐文彬, 张睿, 陆旭兵 2023 物理学报 72 097703 doi: 10.7498/aps.72.20222221 Yuan G L, Wang C H, Tang W B, Zhang R, Lu X B 2023 Acta. Phys. Sin. 72 097703 doi: 10.7498/aps.72.20222221
    [9] Liu K, Liu K, Zhang X, Fang J, Jin F, Wu W, Ma C, Wang L 2024 Chin. Phys. Lett. 41 117701 doi: 10.1088/0256-307X/41/11/117701
    [10] Scott E A, Gaskins J T, King S W, Hopkins P E 2018 APL Mater. 6 058302 doi: 10.1063/1.5021044
    [11] Lu S X, Cebe P 1996 Polymer 37 4857 doi: 10.1016/S0032-3861(96)00376-X
    [12] 何宇亮, 周衡南, 刘湘娜, 程光煦 1990 物理学报 39 1796 doi: 10.3321/j.issn:1000-3290.1990.11.015 He Y L, Zhou H N, Liu X N, Cheng G X, Yu S D 1990 Acta Phys. Sin. 39 1796 doi: 10.3321/j.issn:1000-3290.1990.11.015
    [13] Wang S, Wang C, Li M, Huang L, Ott R, Kramer M, Sordelet D, Ho K 2008 Phys. Rev. B 78 184204 doi: 10.1103/PhysRevB.78.184204
    [14] Wang Y, Fan Z, Qian P, Caro M A, Ala-Nissila T 2023 Phys. Rev. B 107 054303 doi: 10.1103/PhysRevB.107.054303
    [15] 林揆训, 林璇英, 梁厚蕴, 池凌飞, 余楚迎, 黄创君 2002 物理学报 51 863 doi: 10.7498/aps.51.863 Lin K X, Lin X Y, Liang H Y, Chi L F, Yu C Y, Huang C J 2002 Acta. Phys. Sin. 51 863 doi: 10.7498/aps.51.863
    [16] Kittel C 1949 Phys. Rev. 75 972 doi: 10.1103/PhysRev.75.972
    [17] Allen P B, Feldman J L 1989 Phys. Rev. Lett. 62 645 doi: 10.1103/PhysRevLett.62.645
    [18] Zhu X, Shao C 2022 Phys. Rev. B 106 014305 doi: 10.1103/PhysRevB.106.014305
    [19] Qiu R, Zeng Q Y, Han J S, Chen K, Kang D D, Yu X X, Dai J Y 2025 Phys. Rev. B 111 064103 doi: 10.1103/PhysRevB.111.064103
    [20] Luo T, Lloyd J R 2012 Adv. Funct. Mater. 22 2495 doi: 10.1002/adfm.201103048
    [21] Li S H, Yu X X, Bao H, Yang N 2018 J. Phys. Chem. C 122 13140 doi: 10.1021/acs.jpcc.8b02001
    [22] Xi Q, Zhong J, He J, Xu X, Nakayama T, Wang Y, Liu J, Zhou J, Li B 2020 Chin. Phys. Lett. 37 104401 doi: 10.1088/0256-307X/37/10/104401
    [23] Prasai K, Biswas P, Drabold D 2016 Semicond. Sci. Techn. 31 073002 doi: 10.1088/0268-1242/31/7/073002
    [24] 朱美芳 1996 物理学报 45 499 doi: 10.7498/aps.45.499 Zhu M F 1996 Acta Phys. Sin. 45 499 doi: 10.7498/aps.45.499
    [25] Mu X, Wu X, Zhang T, Go D B, Luo T 2014 Sci. Rep. 4 3909 doi: 10.1038/srep03909
    [26] 王学智, 汤雨婷, 车军伟, 令狐佳珺, 侯兆阳 2023 物理学报 72 056101 doi: 10.7498/aps.72.20221581 Wang X Z, Tang Y T, Che J W, Linghu J J, Hou Z Y 2023 Acta Physica Sinica 72 056101 doi: 10.7498/aps.72.20221581
    [27] Yang F H, Zeng Q Y, Chen B, Kang D D, Zhang S, Wu J H, Yu X X, Dai J Y 2022 Chin. Phys. Lett. 39 116301 doi: 10.1088/0256-307X/39/11/116301
    [28] 鲍华 2013 物理学报 62 1 doi: 10.7498/aps.62.186302 Bao H 2013 Acta Phys. Sin. 62 1 doi: 10.7498/aps.62.186302
    [29] Qiu R, Zeng Q Y, Wang H, Kang D D, Yu X X, Dai J Y 2023 Chin. Phys. Lett. 40 116301 doi: 10.1088/0256-307X/40/11/116301
    [30] Isaeva L, Barbalinardo G, Donadio D, Baroni S 2019 Nat. Commun. 10 3853 doi: 10.1038/s41467-019-11572-4
    [31] Simoncelli M, Marzari N, Mauri F 2019 Nat. Phys. 15 809 doi: 10.1038/s41567-019-0520-x
    [32] Harper A F, Iwanowski K, Witt W C, Payne M C, Simoncelli M 2024 Phys. Rev. Mater. 8 043601 doi: 10.1103/PhysRevMaterials.8.043601
    [33] Fiorentino A, Pegolo P, Baroni S 2023 npj Comput. Mater. 9 157 doi: 10.1038/s41524-023-01116-2
    [34] Barbalinardo G, Chen Z, Lundgren N W, Donadio D 2020 J. Appl. Phys. 128 135104 doi: 10.1063/5.0020443
    [35] Fan Z Y, Wang Y Z, Ying P H, Song K K, Wang J J, Wang Y, Zeng Z Z, Xu K, Lindgren E, Rahm J M, Gabourie A J, Liu J H, Dong H K, Wu J Y, Chen Y, Zhong Z, Sun J, Erhart P, Su Y J, Ala-Nissila T 2022 J. Chem. Phys. 157 114801 doi: 10.1063/5.0106617
    [36] Fan Z Y, Zeng Z Z, Zhang C Z, Wang Y Z, Song K K, Dong H K, Chen Y, Ala-Nissila T 2021 Phys. Rev. B 104 104309 doi: 10.1103/PhysRevB.104.104309
    [37] Zhang H, Gu X, Fan Z, Bao H 2023 Phys. Rev. B 108 045422 doi: 10.1103/PhysRevB.108.045422
    [38] Sivaraman G, Krishnamoorthy A N, Baur M, Holm C, Stan M, Csányi G, Benmore C, Vázquez-Mayagoitia Á 2020 npj Comput. Mater. 6 104 doi: 10.1038/s41524-020-00367-7
    [39] Plimpton S 1995 J. Comput. Phys. 117 1 doi: 10.1006/jcph.1995.1039
    [40] Zeng Q Y, Yu X X, Yao Y P, Gao T Y, Chen B, Zhang S, Kang D D, Wang H, Dai J Y 2021 Phys. Rev. Res. 3 033116 doi: 10.1103/PhysRevResearch.3.033116
    [41] Franzblau D 1991 Phys. Rev. B 44 4925 doi: 10.1103/PhysRevB.44.4925
    [42] Le Roux S, Petkov V 2010 J. Appl. Crystallogr. 43 181 doi: 10.1107/S0021889809051929
    [43] Xiang X, Fan H, Zhou Y G 2024 J. Appl. Phys. 135 125102 doi: 10.1063/5.0190047
    [44] Zhang S L, Yi S L, Yang J Y, Liu J, Liu L H 2023 Int. J. Heat Mass Tran. 207 123971 doi: 10.1016/j.ijheatmasstransfer.2023.123971
    [45] Panzer M A, Shandalov M, Rowlette J A, Oshima Y, Chen Y W, McIntyre P C, Goodson K E 2009 IEEE Electron Dev. Lett. 30 1269 doi: 10.1109/LED.2009.2032937
    [46] Lee S M, Cahill D G, Allen T H 1995 Phys. Rev. B 52 253 doi: 10.1103/PhysRevB.52.253
    [47] Chaubey G S, Yao Y, Makongo J P, Sahoo P, Misra D, Poudeu P F, Wiley J B 2012 RSC Adv. 2 9207 doi: 10.1039/c2ra21003g
    [48] Tritt T M 2005 Thermal Conductivity: Theory, Properties, and Applications (Springer Science & Business Media
    [49] Shenogin S, Bodapati A, Keblinski P, McGaughey A J 2009 J. Appl. Phys. 105 034906 doi: 10.1063/1.3073954
    [50] Seyf H R, Henry A 2016 J. Appl. Phys. 120 025101 doi: 10.1063/1.4955420
  • 加载中
图( 7)
计量
  • 文章访问数:  27
  • HTML全文浏览数:  27
  • PDF下载数:  0
  • 施引文献:  0
出版历程
  • 收稿日期:  2025-03-17
  • 刊出日期:  2025-06-05

基于准简谐格林-久保理论结合流体动力学外推方法的非晶二氧化铪导热机制研究

    通讯作者: E-mail: shouhang.li@universite-paris-saclay.fr.; 
    通讯作者: E-mail: xjliu@dhu.edu.cn.
  • 1. 东华大学机械工程学院, 微纳机电系统与集成电路研究所, 上海 201620
  • 2. 巴黎萨克雷大学, 法国国家科学研究中心, 纳米科学与纳米技术中心, 帕莱索 91120, 法国

摘要: 非晶态材料二氧化铪在微电子器件中具有广泛应用, 理解其微观导热机制对于提升电子器件的性能和可靠性至关重要. 以往的研究大多基于分子动力学和单一准简谐格林-久保方法, 难以准确考虑低频振动模式的导热贡献. 本文基于准简谐格林-久保理论, 结合流体动力学外推法, 对不同有序度的非晶二氧化铪结构的热输运机制进行全面研究. 该方法可有效克服单一准简谐格林-久保方法中的有限尺寸问题. 理论预测表明, 非晶二氧化铪的热导率与微观结构有序度呈现弱相关性. 基于模态分析表明, 中低频振动模式对热导率具有显著贡献, 是单一准简谐格林-久保方法低估非晶二氧化铪热导率的主要原因. 同时, 本文基于非谐动态结构因子分离了传播子和扩散子对非晶二氧化铪导热的贡献, 计算表明扩散子在所有非晶二氧化铪结构导热中均占据主导作用. 然而, 传播子的导热贡献仍不可忽略, 其占比可高达20%以上, 且随着有序度的增大而增大.

English Abstract

    • 基于过渡金属氧化物的电阻式随机存储器(RRAM)因具有优良的半导体工艺兼容性、快速开关时间和低能耗等优点[1,2], 是下一代非易失性存储器最有前景的候选器件之一. 非晶二氧化铪(HfO2)具有较强的离子极化和电子极化性能, 在电场作用下可产生较大的电位移, 从而具有较高的介电常数[3,4]. 同时, 因其具有优越的铁电性质[5,6]且与等离子沉积和原子层沉积等工艺技术具有良好的兼容性, 通常被选为制造RRAM器件的开关层材料[79]. 考虑到RRAM器件运行过程中产生的自热效应会严重影响其使用性能和可靠性, 非晶HfO2的微观导热机制引起了人们的广泛关注[10].

      在工业领域, 熔体快速冷却[11]和气相沉积[12]是制造非晶固体电介质材料的两种常用方法. 快速冷却可有效抑制晶体的成核和生长, 在冷却速度足够高时形成非晶态. 在气相沉积法中, 利用化学气相沉积、物理气相沉积、脉冲激光沉积或等离子体增强化学气相沉积分离原子或分子, 然后沉积到低温冷却的基底上形成非晶态. 此方法通过控制冷却速度可以产生不同的非晶形态, 使得不同微观结构的键长、键角、配位数和成环数存在显著差异[13]. 目前在非晶硅中已经发现, 微观结构的有序度对热导率具有明显影响[14,15], 但有序度对非晶HfO2导热的影响仍然未知. 因此, 有必要揭示非晶HfO2的热输运机制, 以实现对其热导率进行有效调控.

      相比晶体材料, 现有关于非晶材料的热传输机制研究相对较少, 理论研究手段较为欠缺. Kittel[16]在20世纪40年代采用传统的声子导热图像唯象解释了非晶材料的温度依赖特性. Allen和Feldman[17](A-F)根据声子形状和特征向量来区分非晶固体中的三类热振动, 即传播子、扩散子和局域子. 他们建立了一个网络模型, 并采用简谐近似来描述非晶材料中的热输运, 解释了在高温下热导率饱和的现象. 然而, A-F模型忽略了非简谐效应, 而非简谐效应会引起频率偏移, 这对于非晶材料热导率的预测非常重要[18,19]. 分子动力学(MD)模拟可以用于构建无序的非晶原子结构, 并考虑任何阶次的非简谐效应, 因此被广泛用于预测非晶聚合物[2022]、半导体[23,24]和氧化物[2527]的热导率. 然而, 由于传统MD方法依赖经验力场, 并忽略了量子效应, 在很大程度上削弱了其准确性. 此外, 采用MD方法获得收敛的热导率计算量巨大, 难以准确揭示低频振动模式对热导率的贡献[28,29].

      近年来发展的线性响应格林-久保方法[30]和Wigner输运方程[31]成功拓展了传统的A-F模型, 为研究非晶材料中复杂的热传输机制提供了统一的理论框架. 这些方法可以同时考虑结构无序性和非谐性的影响, 因此可以在很宽的温度范围内准确预测非晶材料的热导率[32]. 然而, 这些方法的计算复杂度与原子数(N)的关系为$\mathcal{O}\left( {{N^3}} \right)$[33], 仅能计算有限尺寸非晶结构的热导率, 不能充分激发低频长波模式, 因而无法准确考虑低频振动模式的导热贡献. 最近, Fiorentino等[33]发展了流体动力学外推方法, 克服了低频传播模式在有限尺寸模拟中热输运贡献被低估的问题. 在该方法中, 传播子的热导率计算基于连续介质弹性理论, 而扩散子和局域子的热导率则是依据准简谐格林-久保(QHGK)方法[34]. 因此, 采用流体动力外推方法可将有限尺寸系统的热导率外推至无限大尺寸系统所对应的热导率.

      本文采用QHGK方法结合流体动力学外推方法, 全面研究了结构有序度对非晶HfO2晶格热导率的影响. 首先利用MD通过控制淬火速度获得不同有序度的非晶结构, 之后以从MD中提取的力常数为输入, 基于QHGK理论结合流体动力学外推方法计算了其热输运特性, 并讨论了QHGK方法中有限尺寸对热导率的影响. 最后, 分离了不同热振动模式对非晶结构热导率的贡献.

    • 非晶结构通过MD模拟“熔融-淬火-退火”过程获得, 初始原子结构由规则立方体单元沿3个笛卡尔坐标方向分别复制$ n = 3$$ n= 6 $次得到. 所有非晶结构的构建均使用GPUMD软件包[35]进行. 神经演化势函数[36] (NEP)结合了第一性原理计算的高精度和传统势函数的高效性, 本研究中的原子间相互作用由文献[37]开发的NEP机器学习势描述, Zhang等[37]使用了该NEP势函数计算的非晶HfO2的X射线结构因子与实验结果一致, 计算的振动状态密度与密度泛函理论预测相吻合, 因此 该势函数具有较高的可靠性. 未经特殊说明, 本文所有MD模拟的时间步长均为1 fs. MD模拟过程如下: 初始结构首先在20 ps内从300 K升温到4000 K, 之后样品在4000 K的正则系综(NVT)中弛豫50 ps, 然后再以不同的淬火速度降温至300 K. 将淬火所得结构在300 K的等温等压系综(NPT)下进行500 ps退火处理, 选取该系综下弛豫结构的最后一帧得到最终的非晶构型. 上述所有的等温等压模拟均在零压力下进行. 为了避免MD模拟中随机误差对最终结果的影响, 本文采用不同的随机数种子初始化原子速度, 在每个淬火速度下均构建了5个不同的结构. 最终获得的非晶结构的质量密度约为7.69 g/cm3, 与实验测量结果和其他理论预测结果吻合较好[38], 进而侧面验证了本文非晶HfO2结构生成过程的可靠性. 需要注意的是, 本文关注的重点是最终的非晶结构有序度对热导率的影响, 因而暂不讨论中间模拟过程参数(压强、温度等)对非晶结构构建的影响.

    • 本文利用$ \kappa{\mathrm{ALD}}o$软件包[34]结合大规模原子分子并行模拟器(LAMMPS)[39]提取非晶HfO2结构的原子间作用力常数, 之后基于QHGK理论计算非晶结构的热导率. 本文所采用的广义热导率表达式如下[34]:

      其中, $ V $是系统的体积, $ \alpha $$ {\alpha }' $是笛卡尔坐标, $ \mu $$ {\mu }' $表示不同的振动模式. 广义热容$ {C}_{\mu {\mu }'} $可表达为

      式中, $ T $是温度, 玻色-爱因斯坦分布$ {n}_{\mu }= {\left({{\mathrm{e}}}^{\hbar{\omega }_{\mu }/{k}_{{\mathrm{B}}}T}-1\right)}^{-1} $. 广义声子群速度$ {v}_{\mu {\mu }'} $的表达式为

      其中, $ x $$ M $是分别是系统中原子的位置和质量, $ {\varPhi _{i\beta j\beta '}} $=$ \dfrac{{\partial }^{2}{E}_{{\mathrm{P}}}}{\partial {R}_{i\alpha }\partial {R}_{j\alpha }} $是二阶原子间作用力常数, $ {E}_{{\mathrm{P}}} $是势能, $ \boldsymbol{\eta } $是动力学矩阵对应的特征向量. 广义寿命$ {\tau }_{\mu {\mu }'} $是一个类洛伦兹形分布量, 其量化了接近共振的振动模式之间的扩散过程:

      其中, $ {\varGamma }_{\mu } $是振动模态$ \mu $的非简谐线宽, 可以用费米黄金法则计算[34].

    • 将QHGK方法直接应用于原子数大于1000的非晶系统需要消耗较大的计算机内存和较长的计算时间. 为了克服有限尺寸问题, 获得无限大系统的热导率, 本文采用文献[33]提出的流体动力学外推方法. 该方法分离了总热导率中的传播子($ {\kappa }_{{\mathrm{P}}} $)贡献与扩散子和局域子($ {\kappa }_{{\mathrm{D}}} $)贡献:

      将频率低于截止频率$ {\omega }_{{\mathrm{P}}} $的声学平面波视作传播子, 其热导率可表达为

      其中$ b={\mathrm{L}}, {\mathrm{T}} $表示纵向($ {\mathrm{L}} $)和横向(T)声学分支. $ \varGamma $表示线宽. 对于低频长波振动模式, 其线宽满足$ {\varGamma }_{b}\left(\omega \right)={\alpha }_{1}{\omega }^{2}+{\alpha }_{2}{\omega }^{4} $, 其中$ {\alpha }_{1} $$ {\alpha }_{2} $是拟合参数. $ {c}_{b} $是声速, $ {\rho }_{b} = {{\omega }^{2}}/({{\pi }^{2}{c}_{b}^{3}}) $是态密度. 扩散子的热导率对体系尺寸不敏感, 采用(7)式进行计算:

      式中, Θ是Heaviside-theta函数, 其用于区分传播模式和非传播模式的贡献. 传播子、扩散子和局域子可根据其振动特征加以区分, 这些特征可以由非谐动态结构因子来反映. 简谐系统的动态结构因子[40]表达式可通过引入非简谐效应加以推广:

      式中, $ \langle\mu |\boldsymbol{Q}, b\rangle $表示第$ \mu $个简正模在由波矢$ \boldsymbol{Q} $所表征的振动方向上的投影.

      本文采用LAMMPS, 采用有限差分法分别提取了324个原子非晶HfO2结构的二阶和三阶原子间作用力常数以及2592个原子的非晶HfO2结构的二阶原子间作用力常数. 在提取力常数过程中, 原子位移设置为1×10–6 Å. 随后将原子结构和相应的力常数作为QHGK求解器$ \kappa{\mathrm{ALD}}o $的输入以求解热导率. 为了计算(6) 式中的$ {\kappa }_{{\mathrm{P}}} $值, 需要拟合出$ {\varGamma }_{b}\left(\omega \right) $的函数表达式. 首先将324和2592个原子的非晶结构计算出的非谐动态结构因子((8)式)拟合成下面的洛伦兹函数形式:

      进而获得声速$ {c}_{b} $和线宽$ {\varGamma }_{b} (Q) $. 随后采用线性色散, 即$ {\omega }_{b} = {c}_{b}Q $, 得到线宽与频率的函数关系$ {\varGamma }_{b} (\omega) $, 将其拟合成四次多项式即可获得拟合系数$ {\alpha }_{1} $$ {\alpha }_{2} $的数值.

    • 图1可以看出, 在20 ps的熔融过程中, 体系的势能迅速上升, 所对应的均方位移也逐渐增大; 随后在4000 K下弛豫50 ps, 此时体系的势能在–3050 — –3000 eV之间上下波动并且这段过程中的均方位移的相对值跨度很大, 这表明体系内的原子位置变化剧烈; 之后体系以1×1012 K/s的淬火速度降温至300 K, 体系的势能逐渐下降而均方位移变化放缓, 这表明体系内的原子位置变化逐渐变小; 最后在300 K的NPT系综下退火500 ps, 体系的势能在–3235.6 — –3238.3 eV之间产生很小的波动, 均方位移的累计值变化范围为2.4574—2.4710 Å2, 这说明淬火后的结构在退火500 ps后已经足够稳定.

    • 考虑到非晶态结构不具备长程有序特征, 本文采用短程和中程有序来量化非晶结构的有序度. 短程有序可由径向分布函数(RDF)来描述, 如图2(a)所示. 非晶HfO2的第1个RDF峰位于2.13 Å处, 表征了Hf-O原子对的平均键长, 该数值与文献[38]中保持一致. 第1个RDF峰明显大于其他峰, 这是因为非晶HfO2的长程有序性相对较弱. 如图2(a)插图所示, 所有非晶结构的RDF的大小随着淬火速度的减小而增大, 这表明短程有序度随着淬火速度的减小而增大. 非晶体系的中程有序可以用原子成环分布[41]来衡量, 图2(b)所示的环分布采用最短路径算法[42]计算得出. 在晶体HfO2中, 八元环因其具有较高的能量稳定性而占主导地位, 但在非晶结构中其比例相对较小且随淬火速度的增大而明显下降. 非晶结构中六元环和十元环的比例随着淬火速度的增大而逐渐增大, 这表明中程有序性逐渐变差.

    • 图3所示为不同淬火速度下采用单一QHGK方法和QHGK方法结合流体动力学外推预测的非晶HfO2热导率. 虽然单一QHGK方法能够反映出不同无序结构的热导率差异, 但是该方法无论在低淬火速度5×1010 K/s还是高淬火速度5×1012 K/s下均明显低估了热导率. 这主要是因为单一QHGK方法只能考虑有限尺寸大小的体系, 低估了低频长波振动模式对导热的贡献, 而流体动力学外推方法可以克服有限尺寸的问题. 从图3可以得出, 无论是采用单一QHGK还是QHGK结合流体动力学外推法, 它们预测的热导率均随着温度的增大而逐渐饱和. 此外, 两种方法预测的热导率的差异逐渐变小, 这主要归因于离域模式的平均自由程随温度增大迅速下降, 有限模拟尺寸对单一QHGK方法预测结果的影响降低. 采用单一QHGK方法预测的有序度对热导率的影响在高温下几乎消失, 而QHGK结合流体动力学外推法计算的热导率表明有序度在高温下对热导率仍存在一定影响. 这反映出即使在高温下仍需采用流体动力学外推法预测热导率.

      图4所示为非晶HfO2随温度变化的热导率. 热导率随温度升高而增大, 这主要归因于晶格热容的增大. 随着温度的进一步升高, 由于热容的增大放缓且热载流子弛豫时间明显下降, 热导率增幅明显降低. 文献[10,43,44]报道的晶体HfO2的室温热导率为2.60—11.43 W/(m·K), 明显高于本文预测的非晶结构的热导率. 整体上, 非晶HfO2热导率随着淬火速度的增大有所下降, 然而, 相较非晶硅而言非常的微弱[14]. 不同研究小组测得的HfO2热导率实验值存在较大差异[10,4547]. 但总的来说, 这些数值都分散在本文的理论预测结果附近, 这进一步验证了本文预测方法的准确性. 图4中实验数据的差异可能归因于所采用的非晶样品微观结构以及实验测量方法的不同.

      不同非晶结构的热输运性质与微观输运特性密切相关, 包括每种振动模式的弛豫时间、模态扩散率和热容, 如图5所示. 如图5(a)所示, 非晶HfO2不同振动模式对热容的贡献随着频率的增大而明显下降. 由于同一材料的不同非晶结构的原子组成保持一致, 因此它们的热容几乎是重叠的. 不同非晶HfO2结构的振动模式弛豫时间见图5(b). 弛豫时间随着淬火速度的降低而逐渐增大, 这表明在淬火速度较大的无序结构中, 非谐性更强, 热载流子散射更多. 然而不同结构的非晶HfO2的弛豫时间的差异并不那么显著, 这就导致了图4所示的热导率与淬火速度呈现出弱相关性. 在QHGK方法中, 模态扩散率可表示为$ {D_\mu } = \displaystyle\frac{1}{3}\sum\nolimits_{^{\mu '}} {{D_{\mu \mu '}}} = \displaystyle\frac{1}{3}\sum\nolimits_{^{\mu '}} {|{\nu _{\mu \mu '}}{|^2}{\tau _{\mu \mu '}}} $. 不同非晶HfO2结构的模态扩散率如图5(c)所示, 0—5 THz振动模式相对应的模态扩散系数明显大于其他振动模式. 这种现象与晶体结构中的情况不同, 在晶体结构中, 低频模态是主要的热传导因素[48], 中频和高频模态的贡献可以忽略不计. 由于非晶HfO2在一定程度上保留了短程有序, 因而还能部分保留晶体的高频振动特征, 振动模态密度在20 THz附近的峰可视为晶体HfO2光学声子模态的残留. 而模态扩散率与振动模态密度相关, 所以在20 THz处有一个明显的峰值. 特别地, 非晶HfO2中20 THz以上振动模式的模态扩散率呈现明显下降趋势, 这表明这部分振动模式属于局域子[49]. 模态参与比$ p(\sigma) = { \Big(N\displaystyle\sum\nolimits_{n = 1}^N {{{\left| {{{\text{e}}^n}(\sigma )} \right|}^4}} \Big)^{ - 1}} $可以定量表征振动模式的空间局域化程度[50], 其中$ {{\mathrm{e}}}^{n} $是由作用于原子$ \sigma $的第$ n $个特征向量的3个笛卡尔坐标分量组成的向量. 如图5(d)所示, 非晶HfO2中的振动模式主要由离域模式组成(参与比p ~ 0.3), 而当频率高于20 THz时, 参与比明显下降, 这表明振动模式已转变为局域子. 随着淬火速度的降低, 低频振动模式参与比有一定程度增大, 中频模式参与比略微升高, 而高频振动模式的参与比几乎不发生变化. 结合图5(a), (c)可以看出, 非晶HfO2导热主要由中低频振动模式贡献, 这是图3中单一QHGK方法明显低估热导率数值的主要原因.

      不同振动模式的热传输行为与其振动特性密切相关, 根据非谐动态结构因子, 频率低于截止频率$ {\omega }_{{\mathrm{P}}} $的振动模式被归类为类声子输运(传播子), 而频率高于截止频率的模式则被视为非传播模式(扩散子和局域子). 如图6所示, 非晶HfO2的非谐动态结构因子在低频区呈线性变化, 故本文将截止频率$ {\omega }_{{\mathrm{P}}} $= 3 THz以下的振动模式归类为传播子. 这里只区分振动模式中的传播子和扩散子对导热的贡献, 因为上文已分析非晶HfO2中局域子对热导率的贡献可以忽略不计.

      依据上述区分热载流子的方法, 本文分离了不同淬火速度下非晶HfO2中传播子和扩散子对热导率的贡献. 如图7所示, 对于所有的非晶HfO2结构, 传播子对热导率的贡献均在20%以上, 明显高于Shenogin等[49]估计的非晶氧化物中传播子对导热的贡献. 这进一步说明了非晶结构热导率预测中考虑有限尺寸问题的重要性. 随着淬火速度的降低, 传播子对热导率的贡献比例逐渐增大. 这是因为随着淬火速度的降低, 非晶体系中的微观结构变得更加有序(如图2所示), 从而导致更多的热载流子向类声子输运模式转变(如图5(d)所示), 即低频区域的振动模式在热输运过程中的贡献越来越高. 值得注意的是, 在所有的非晶HfO2结构中, 扩散子对导热贡献均占据主导地位(高于70%).

    • 本文采用准简谐格林-久保理论并结合流体动力学外推方法, 系统研究了微观原子结构对非晶HfO2晶格热导率的影响. 该方法可有效克服单一准简谐格林-久保方法研究非晶导热存在的有限尺寸问题. 利用分子动力学模拟, 通过控制淬火速度生成了不同有序度的非晶结构, 理论预测的非晶HfO2热导率数值能在较宽的温度范围内匹配实验数据. 与非晶硅不同, 非晶HfO2的热导率与微观结构有序度呈现出弱相关性. 基于模态分析表明, 中低频振动模态对热导率具有显著贡献, 这是单一准简谐格林-久保方法低估非晶HfO2热导率的主要原因. 本文基于非谐动态结构因子分离了传播子和扩散子对非晶HfO2的导热贡献, 结果表明传播子的贡献不可忽略, 且其重要性随着结构有序度的增大而逐渐增大.

    参考文献 (50)

目录

/

返回文章
返回