基于声子输运理解非平衡分子动力学

分子动力学是微纳尺度传热研究中常用的计算方法。因为其使用简单、容易展示，在微纳尺度传热研究中起着至关重要的作用。分子动力学计算热导率的主流方法有两种，一种是平衡态分子动力学：运用涨落耗散的Green-Kubo公式计算，后面简称EMD。第二种是非平衡分子动力学方法：通过引入热源和冷源，在体系中形成温度梯度来计算热导率，后面简称NEMD。这两种方法各有优劣，但是一个核心问题在于，其背后的声子图像一直以来是不甚清晰的。

 

本文主要讨论NEMD。众所周知，如何建立非平衡对于计算结果影响极大。其中的主要问题包括：

1．体系长度、热源长度、热源类型很大程度上会影响计算结果，导致了对于同样的原子模型，同样的相互作用势，算出来的热导率不一样。

2. 提取热导率的方式也会影响计算结果。有时在热源附近会有明显的非线性现象，因此一般用线性区间的温差和温度梯度来提取热导率。因为计算有尺寸效应，因此一般使用线性外推的方式来得到体块材料的热导率。

3. NEMD两侧边界的形貌会明显影响界面热阻的计算结果(Liang and Keblinski Phys. Rev. B 90, 075411,2014)。

 

上述问题直接导致了对NEMD计算结果缺乏真正的理解。那么，NEMD是否有统一的基于声子图像的物理解释呢？

 

最近，上海交通大学密西根学院鲍华课题组与美国橡树岭国家实验室冯天力博士，芬兰阿尔托大学樊哲勇博士，上海交通大学机动学院顾骁坤副教授以及美国普渡大学的Mark Lundstrom教授等人合作，在Phys. Rev. B上发表了一篇题为“Unification of Non-equilibrium Molecular Dynamics and Mode-resolvedPhonon Boltzmann Transport Equation for Thermal Transport Simulations“的文章。该研究通过定量的比较了NEMD和非灰（nongray）声子玻尔兹曼输运方程（简称BTE）两种方法，明确了NEMD的声子图像。

该论文以Tersoff作用势下的硅作为研究对象。为保证两个模拟完全可比，在BTE的声子弛豫时间、群速度等输入参数也完全从该相互作用势提取。在BTE计算中，采用了经典的声子分布函数。在相同的体系下，我们发现，NEMD和BTE的模拟结果在数值上完全一致（如图1）。实现上述结果的条件是：确保两种方法的热源加载方式相同，即Noose Hoover热源对应BTE中的均匀内热源，而Langevin热源则对应于BTE中的定温边界条件。

图1.(a) Langevin热源的NEMD和定温边界条件的BTE的温度与热流分布对比。(b) NHC热源的 NEMD和均匀内热源的BTE的温度与热流分布对比。(c) Langevin热源的 NEMD和定温边界条件的BTE的声子支温度对比。(d) NHC热源的 NEMD和均匀内热源的BTE的声子支温度分布对比。

 

基于该对比，发现NEMD中Langevin这种随机性控制的热源和Nose Hoover这类确定性控制的热源，在声子图像上是完全不同的。Langevin热源对应于BTE中的定温边界条件，因此热源和样品的界面处会发射热源温度的声子，而所有的声子在进入热源后，都会迅速达到平衡，而且在热源内部声子也不会与边界发生散射。而对于Nose Hoover热源以及Velocity rescaling热源，对应于BTE中均匀加热的内热源，其本质上没有控制热源中的声子的温度，只能通过调整内热源的能量大小来维持热源的平均温度，因此声子进入其中，不会迅速达到热平衡，而会在热源内部继续输运，有些声子还会与边界散射。（如图2）

图2. (a) 采用Langevin热源的NEMD模拟的物理图像。(b) 采用Nose Hoover热源或Velocity rescaling热源的 NEMD模拟的物理图像

 

基于这些理解，前言中提到的NEMD的问题也就迎刃而解了。我们认为，

1. 体系的长度会造成热导率结果的不同，就是声子的弹道输运导致的。在本公众号早期的一篇文章中有介绍过，这个效应也在NEMD和BTE中都有过很多研究。
2. 不同的热源类型会造成不同的结果，这是因为热源加载方式本质上是不同的。Langevin热源保证了内部的声子都处于平衡态，而Nose Hoover、Velocity rescaling等热源内部的声子都偏离了平衡态（如图1(c)和1(d)）。
3. 如何提取热导率？基于上述理解，就可以得到一个比较推荐的提取方案：可以采用Langevin热浴建立非平衡，然后利用热浴的温差（而非温度梯度）除以热流的方式来提取热导率。这样的数值即为名义热导率（Apparent thermal conductivity）。这种定义在BTE和朗道方法中一直被广泛采取。这种方法计算出的结果在弹道极限下可以和Atomistic Green’s function的结果一致（Li et al, Journal of Chemical Physics, 151, 234105, 2019）。而Nose Hoover中，无法明确的定义温差，而且计算结果和热浴长度有关。
4. 对于边界形貌这个问题，采用Langevin热源，声子不会与边界发生散射，因此NEMD两侧边界的形貌也不会对结果造成影响；采用Nose Hoover热源或者Velocity rescaling热源，声子与边界发生散射，NEMD两侧边界的形貌均对结果造成影响（如图3(c)和3(d)）。

图3. (a)不同Langevin热源长度的NEMD模拟温度场结果。(b)不同Nose Hoover热源长度的NEMD模拟温度场结果。(c)采用Langevin热源不同边界形貌的NEMD模拟温度场结果。(d)采用Nose Hoover热源不同边界形貌的NEMD模拟温度场结果。(e)不同边界形貌的NEMD模拟域，和图(c) 和(d)中的图例对应。

 

论文信息：

Yue Hu, Tianli Feng, Xiaokun Gu, Zheyong Fan, Xufeng Wang, Mark Lundstrom, Som S. Shrestha, Hua Bao*.Unification of nonequilibrium molecular dynamics and the mode-resolved phonon Boltzmann equation for thermal transport simulations. Phys. Rev. B 101, 155308,2020