揭秘CuCl超低热导率：四声子散射与温度重正化的关键作用

1. 项目概述：为何要深挖CuCl的热导率？

在材料科学和凝聚态物理的交叉领域，热输运性质的研究从来都不是一个孤立的课题。它直接关系到热电材料的转换效率、电子器件的散热能力，以及热障涂层的服役寿命。传统上，我们倾向于认为结构简单、原子质量差异大的二元半导体，比如硅化锗（GeSi）或砷化镓（GaAs），通常具有较高的热导率。这是因为大的质量差会导致声子谱中光学支和声学支分离，减少了声子间的散射通道，热量传导更顺畅。然而，氯化亚铜（CuCl）这个材料却是个“异类”。它拥有极其简单的闪锌矿结构，铜和氯的原子质量也相差不小，但实验测得其室温热导率却低得惊人，仅有约0.7到0.85 W m⁻¹ K⁻¹，比许多结构复杂的玻璃还要低。这个“简单结构，超低热导”的矛盾现象，就像一个物理谜题，吸引了众多研究者的目光。

早期的理论计算试图解开这个谜团，但结果却莫衷一是。有些研究严重高估了热导率，有些则低估了。这种分歧指向了一个核心问题：对于像CuCl这样具有强非谐性的材料，传统的计算框架可能“力不从心”。非谐性，简单说就是原子间的相互作用力不严格遵循胡克定律（力与位移成正比），当原子振动偏离平衡位置较大时，势能曲线不再是完美的抛物线。这种非谐性会引发强烈的声子-声子散射，是降低热导率的关键。但传统计算通常只考虑三声子散射过程，并且使用绝对零度（0 K）下的原子间力常数，这无疑忽略了两大关键因素：一是更高阶的四声子散射，二是温度升高导致的原子振动加剧对力常数本身的修正，即力常数重正化。

因此，我们这项工作的出发点非常明确：构建一个更完备、更精确的理论框架，来“复盘”CuCl的热输运行为。我们不仅要解释它为何如此之“冷”（低热导），还要探究一个更反直觉的现象：通常情况下，对材料施加压力会压缩晶格，增强原子间相互作用，使声子“变硬”、传播更快，从而提升热导率。但实验发现，CuCl的热导率在相变发生前，竟然随着压力增加而单调下降。这背后的物理机制是什么？是哪些微观参量在“唱反调”？为了回答这些问题，我们采用了“组合拳”策略：用第一性原理计算提供高精度的基础物理图像，用温度依赖有效势（TDEP）方法严格处理非谐性与温度效应，并引入机器学习力场（MLFF）来攻克高压下计算量巨大的难题。最终，我们不仅成功复现了实验值，还像侦探一样，从声子谱、散射率、群速度等微观量中，找到了导致热导率异常行为的“元凶”。

2. 方法论全景：从第一性原理到热导率计算的完整链条

要准确计算晶格热导率，尤其是对于强非谐性材料，不能简单地套用公式。我们需要一个自洽、严谨且计算上可行的流程。下图概括了我们从原子尺度出发，最终得到宏观热导率的核心工作流：

flowchart TD A[第一性原理计算<br>（VASP, PBEsol泛函）] --> B[原子运动轨迹与受力] B --> C{温度依赖有效势<br>（TDEP）方法} C --> D[提取重正化的<br>谐波与非谐波力常数] D --> E[构建声子谱与散射矩阵] E --> F[求解玻尔兹曼输运方程<br>（包含三/四声子过程）] G[机器学习力场<br>（MLFF）训练与验证] --> H[高效生成高压下的<br>原子受力数据] H --> C F --> I[获得晶格热导率 κl<br>及其对角与非对角贡献] C --> J[微观机制分析<br>（群速度、散射率、态密度等）] I --> K[与实验数据对比验证] J --> K

这个工作流的核心在于几个关键方法的选择与耦合，下面我们来逐一拆解。

2.1 理论基础：超越三声子的热输运方程

晶格热导率（κl）的计算核心是求解声子玻尔兹曼输运方程（BTE）。对于强非谐性材料，必须超越传统的单弛豫时间近似和仅包含三声子的微扰论。我们采用了基于模式耦合理论的非微扰框架，该框架通过声子关联函数来捕捉对热输运有贡献的集体和相干效应。

总热导率张量 κ 由对角贡献（κd）和非对角贡献（κnd）两部分组成：κ = κd + κnd。

• 对角贡献（κd）：这部分来源于声子模式自身的输运，其表达式包含了所有声子模式的群速度、热容以及关键的散射矩阵 Σ的逆。散射矩阵的对角元对应着声子模式的散射率 Γλ，它包含了三声子、四声子以及同位素散射的贡献。计算 κd 的挑战在于需要对角化这个巨大的散射矩阵，我们采用了迭代求解法来高效处理。
• 非对角贡献（κnd）：这部分描述了不同声子模式之间的耦合对热流的贡献，即相干输运。对于具有复杂晶格、许多原子原胞和密集声子分支的材料（如某些笼状化合物或复杂氧化物），这部分贡献可能非常显著。其表达式涉及不同模式间的非对角群速度和耦合散射率。

在CuCl这个案例中，由于其简单的闪锌矿结构（每个原胞仅2个原子），声子分支少且间隔大，模式间的耦合很弱。我们的计算证实，其非对角贡献在室温下仅为约0.04-0.06 W m⁻¹ K⁻¹，相对于总热导率可以忽略不计。这意味着对于CuCl，热量主要由独立的声子模式携带，相干效应很弱。这简化了分析，让我们可以更聚焦于对角贡献中散射过程的物理。

注意：方法选择的关键。对于简单晶体，忽略非对角贡献是合理的近似，能大幅节省计算资源。但对于复杂材料（如SKU、Clathrate等），非对角贡献可能占总热导的相当比例，此时必须采用完整的模式耦合理论或Wigner公式进行计算，否则会严重低估热导率。

2.2 核心输入：如何获得“温度感知”的原子间力常数？

无论是计算对角还是非对角贡献，都需要最基础的输入：谐波和非谐波原子间力常数（IFCs）。IFCs本质上描述了原子偏离平衡位置时，感受到的恢复力的大小，是连接原子微观运动和宏观声子性质（如频率、群速度）的桥梁。

对于强非谐性材料，最大的陷阱在于使用静态的、0 K下的IFCs。温度升高时，原子在平衡位置附近振动的振幅增大，其平均所处的势能面区域与绝对零度时不同。这会导致有效的“力常数”发生改变，即温度诱导的IFC重正化。忽略这一点，就如同用冰点的水分子相互作用来预测沸点水的性质，必然产生巨大偏差。

我们采用TDEP（Temperature Dependent Effective Potential）方法来获取温度依赖的IFCs。其核心思想是“用有限温度下的原子运动来反推有效势”。具体步骤是：

1. 生成原子构型：在目标温度下，通过分子动力学或随机采样，生成一系列包含原子热扰动的超胞构型。
2. 第一性原理计算力：对每一个构型，使用密度泛函理论（DFT）计算每个原子所受的力。
3. 迭代拟合：将原子位移和对应的受力数据，拟合到一个有效势能模型中（通常展开到位移的二阶、三阶、四阶项）。TDEP通过迭代优化，使得该有效势给出的力与DFT计算的力最接近，从而得到该温度下“重正化”后的各阶IFCs。

在我们的计算中，我们构建了包含250个原子（5×5×5超胞）的模型来提取IFCs。对于三阶和四阶非谐力常数，分别设置了6 Å和4 Å的��断半径，在计算精度和效率间取得了平衡。所有DFT计算均使用VASP软件包完成，并经过仔细的收敛性测试。

2.3 高压研究的利器：机器学习力场（MLFF）

研究压力对热导率的影响，意味着需要在多个压力点重复上述流程：每个压力点都需要进行第一性原理分子动力学采样和大量的力计算，以得到该压力下温度依赖的IFCs。这对计算资源是巨大的挑战。

为此，我们引入了VASP的在线机器学习力场（on-the-fly MLFF）。其工作流程是：

1. 初始采样与训练：在少数几个压力-温度条件下，运行短时间的第一性原理分子动力学，同时收集原子位置和受力的数据对。
2. 力场构建：利用这些数据，实时训练一个高维神经网络势函数（或类似模型），使其能够以接近DFT的精度，预测原子受力。
3. 验证与生产：对训练好的MLFF进行严格验证，包括比较能量、力、声子谱乃至最终的热导率与DFT基准结果的差异（如图5、图6所示）。验证通过后，便可用这个MLFF替代昂贵的DFT，快速生成高压下大量构型的受力数据，供TDEP提取IFCs。

我们验证的结果非常令人鼓舞：MLFF预测的能量和力与DFT结果相比，均方根误差（RMSE）极低（能量~0.0004 eV/atom，力~0.0065 eV/Å）。更重要的是，由MLFF IFCs计算出的300 K热导率，与DFT直接计算的结果偏差在6%以内。这证明了MLFF不仅能复现原子间作用，更能准确传递到宏观热物性，是进行高压高通量研究的可靠工具。

3. 计算细节与参数选择：魔鬼在细节中

理论框架搭建好后，具体的计算参数选择直接决定了结果的可靠度。这里分享一些关键的选择和背后的考量。

3.1 交换关联泛函之争：PBEsol的胜出

DFT计算中，交换关联泛函的选择对晶格常数、声子谱等基础性质有决定性影响。我们对比了LDA、PBE和PBEsol三种泛函：

• PBE：给出了最接近实验值的晶格常数（5.43 Å vs. 实验 5.424 Å）。
• PBEsol：晶格常数略小（5.30 Å），但其计算的声子谱，尤其是在对热导率贡献最大的声学支区域，与低温实验数据符合得最好。
• LDA：严重低估晶格常数（5.22 Å），声子谱偏差较大。

这个对比揭示了一个重要现象：最准的晶格常数不等于最准的声子谱和热导率。对于CuCl，PBEsol泛函更好地描述了其电子结构和由此衍生的原子间力，因此在后续所有热导率计算中，我们均采用PBEsol。作为验证，我们也用PBE和LDA计算了热导率，发现它们分别高估了约500%和300%，这进一步确认了PBEsol的适用性。

3.2 收敛性测试：确保结果稳健

第一性原理计算中，收敛性测试是避免“垃圾进，垃圾出”的保险栓。我们进行了系统测试：

• 平面波截断能：设置为500 eV，确保总能量收敛。
• k点网格：对于超胞力计算，采用Γ点采样。对于电子结构计算，使用了更密的网格。
• 超胞尺寸与IFC截断：5×5×5的超胞（250原子）结合6 Å/4 Å的截断半径，足以使三阶、四阶IFCs收敛。我们在补充材料中展示了热导率随截断半径和q点网格密度的收敛情况。
• TDEP采样：确保生成的原子位移构型足够多，以覆盖相空间，从而稳定地拟合出高阶IFCs。

实操心得：IFC截断的艺术。设置非谐IFC的截断半径是个经验活。截断太小，会丢失重要的长程相互作用；截断太大，不仅计算量激增，拟合过程还可能因参数过多而不稳定。我们的策略是：先从较小的截断开始计算热导率，逐步增大截断，直到热导率的变化在可接受的误差范围内（例如<2%）。对于CuCl，4-6 Å是一个较好的平衡点。

3.3 四声子散射的纳入：从三体到四体相互作用

传统计算只考虑到三声子过程（一个声子衰变成两个，或两个声子合并成一个）。但对于强非谐性材料，四声子散射（涉及四个声子的碰撞过程）可能变得非常重要。在TDEP框架中，一旦提取了四阶IFCs，就可以直接计算四声子散射率。

计算四声子散射的挑战在于其相空间（可能的散射通道数）远大于三声子过程，直接计算非常耗时。我们采用了基于选择规则和能量-动量守恒的近似方法，高效地筛选出对热阻贡献最大的四声子散射通道。结果显示，对于CuCl，四声子散射的强度比三声子散射高出几个数量级，这是导致其超低热导率的最直接原因。

4. 结果分析与物理洞察：揭开CuCl热导率之谜

基于上述严谨的计算框架，我们得到了一系列与实验高度吻合的结果，并从中提取了深刻的物理洞察。

4.1 声子谱：温度与压力的烙印

声子谱是理解热输运的起点。我们的计算显示：

• 温度效应：随着温度从0 K升至300 K，CuCl的声子模式整体发生“硬化”（频率升高），尤其在Γ点和X点附近的高对称点处最为明显。这直接体现了温度重正化对IFCs的影响。相比之下，晶格热膨胀对声子谱的影响微乎其微。
• 压力效应：施加压力后，大多数声子模式如预期般硬化。然而，一个反常现象出现了：在X点和L点，横向声学（TA）模式竟然发生了软化（频率降低）。这通常意味着在这些特定的晶格振动模式下，原子间的有效相互作用反而变弱了。这种软化与CuCl中特殊的、反键合主导的电子结构有关（后文详述），是理解其热导率随压力下降的关键线索之一。

4.2 热导率：重正化与四声子散射的压倒性影响

热导率计算的结果清晰地揭示了两个核心因素的巨大作用：

1. IFC重正化：如果固执地使用0 K的IFCs来计算300 K的热导率，即使考虑四声子散射，结果（0.27 W m⁻¹ K⁻¹）也远低于实验值。而使用300 K重正化后的IFCs，计算结果升至0.77 W m⁻¹ K⁻¹，与实验值（0.75-0.85）高度一致。这说明，温度对原子间“弹簧常数”的软化效应，是准确预测热导率的前提。
2. 四声子散射：即使使用了正确的300 K IFCs，如果只考虑三声子散射，热导率会被严重高估（4.75 W m⁻¹ K⁻¹）。一旦加入四声子散射，热导率骤降84%，降至与实验相符的水平。这定量地证明了四声子散射在CuCl中不是可忽略的修正，而是主导热阻的机制。

温度依赖关系也很有趣：仅有三声子散射时，κl ∝ T^{-0.55}，偏离了典型的T^{-1}关系。纳入四声子散射后，温度依赖加强至κl ∝ T^{-1.16}，更接近常见绝缘体的行为。这表明四声子散射不仅降低了热导率的绝对值，还改变了其随温度变化的“斜率”。

4.3 模式分解：谁在主导热阻？

通过计算每个声子模式对热导率的贡献（谱热导率），我们发现：

• 主导频段：低于1.5 THz的低频声子贡献了绝大部分热导。这部分主要是声学支声子。
• TA声子的角色转变：仅考虑三声子散射时，TA模式贡献了约51%的热导。但当引入四声子散射后，TA的贡献骤降至36%。这是因为TA声子，尤其是那些在X和L点附近平坦、色散小的模式，其四声子散射率异常之高（比三声子散射高几个量级），被强烈地“扼杀”了其输���能力。
• LA声子的意外贡献：纵向声学（LA）模式贡献了约50%的热导，其中20%来自高于1.5 THz的频率。这是因为在三声子和四声子散射率存在竞争的区域（~2-3 THz），净散射率反而出现了一个相对低谷，使得这些LA声子得以较长程地���播热量。

这个微观图像告诉我们，CuCl的低热导率是TA声子被四声子过程“重点打击”，以及LA声子在特定频段“侥幸逃生”共同作用的结果。

4.4 压力下的反常行为：群速度与散射率的博弈

通常，压力使晶格紧缩，声子速度加快，热导率上升。但CuCl的热导率随压力单调下降。我们的计算完美复现了这一趋势（图7b），并借助MLFF高效揭示了其机理：

1. TA模式群速度下降：压力下，反常软化的TA模式，其群速度（vg）显著降低。群速度是声子输运热量的“车速”，车速慢了，热导率自然下降。
2. 四声子散射率增强：压力同时显著增强了低频TA声子的四声子散射率。这意味着声子“撞车”（散射）更频繁了，平均自由程缩短。
3. 谐波与非谐效应的拆解：我们做了个“控制变量”计算：固定谐波（二阶）IFCs在0 GPa的值，只让非谐（三、四阶）IFCs随压力变化；以及反过来，固定非谐IFCs，只让谐波IFCs变化。结果发现，固定谐波IFCs时，热导率下降得更快（图8a）。这表明，压力下热导率下降的主要驱动力，来自于谐波性质的变化（即声子谱软化导致的群速度下降），而非非谐性（散射）的进一步增强。当然，两者共同作用才完整再现了实验趋势。

4.5 电子结构根源：弱键与反键合态

一切宏观性质的根源都在于电子结构。我们计算了CuCl的电子态密度（DOS）和晶体轨道哈密顿布居（COHP）。分析发现，在费米能级附近，Cu的d轨道和Cl的p轨道存在微弱的杂化，但主要形成的是反键合态。反键合态占据意味着原子间的成键较弱。计算得到的Cu原子较大的均方位移也证实了这一点。这种弱的、离子性主导的键合，是材料强非谐性的电子层次起源。

施加压力后，这些反键合态向低能方向延伸，积分反键合态密度（ICOHP_AB）增加，表明键合进一步减弱。这从电子结构上解释了为何TA模式会在压力下软化：剪切振动更容易发生在键合较弱的体系中。

5. 常见问题、挑战与解决方案实录

在实际操作这套复杂计算流程时，会遇到各种“坑”。这里记录一些典型问题和我们的应对策略。

5.1 计算不稳定与发散问题

• 问题1：TDEP拟合高阶IFCs时发散或不稳定。

  • 原因：采样构型不足或构型代表性不够；位移幅度太大，超出了势能展开的有效范围；截断半径设置不当，导致拟合参数过多。
  • 解决方案：
    1. 增加采样构型数量，确保覆盖平衡位置附近的玻尔兹曼分布。
    2. 控制原子位移的均方根，使其与目标温度下的热振动幅度匹配。通常可以通过调整分子动力学的初始温度或使用更精细的采样方案来实现。
    3. 谨慎选择非谐IFC的截断半径。可以先从较小的值开始，观察热导率是否收敛，而不是盲目追求大截断。
    4. 在TDEP拟合中启用正则化（regularization）选项，防止过拟合。

• 问题2：计算四声子散射时内存或时间消耗爆炸。

  • 原因：四声子散射的相空间随体系原子数和截断半径呈指数增长。
  • 解决方案：
    1. 利用晶格的对称性大幅减少需要独立计算的q点。
    2. 采用能量和动量守恒的严格选择规则，预先过滤掉不可能发生的散射过程。
    3. 对于大型体系或高精度计算，考虑使用近似方法，如只计算对热阻贡献最大的低频声子区域的四声子过程，或采用基于群论和相空间积分的估计方法。

5.2 MLFF训练与验证陷阱

• 问题：MLFF力场在训练集上表现完美，但在预测新构型（特别是高压、高温）时误差很大。
  • 原因：训练数据没有充分覆盖目标压力-温度相空间；存在外推。
  • 解决方案：
    1. 主动学习：采用“on-the-fly”方式，当MLFF对某个新构型的预测不确定性（如模型方差）超过阈值时，自动调用DFT计算该构型，并将其加入训练集，更新力场。
    2. 分层采样：在构建训练集时，不仅要在环境压力下采样不同温度，还要在目标压力范围内采样多个压力点。确保训练数据在相空间中是均匀分布的。
    3. 多指标验证：不要只满足于能量和力的误差小。必须验证声子谱是否与DFT一致，这是力场能否描述动力学性质的关键。更进一步，像我们做的一样，用MLFF计算一个宏观可观测物理量（如热导率），与DFT基准对比，这是最严格的验证。

5.3 结果分析与物理诠释误区

• 问题：看到热导率计算值与实验吻合，就认为模型完全正确。
  • 警示：数值上的吻合可能是多种误差抵消的结果。必须进行交叉验证和微观机制分析。
  • 正确做法：
    1. 参数敏感性测试：改变截断能、k点密度、超胞大小、IFC截断等，观察热导率的变化是否在可接受范围。
    2. 方法对比：如果可能，用另一种独立的方法（如平衡分子动力学EMD或非平衡分子动力学NEMD）进行交叉验证。
    3. 分解贡献：像我们一样，分别计算三声子和四声子散射的贡献、对角和非对角贡献、不同声子支的贡献。这不仅能验证模型内部是否自洽，更能深入理解物理机制。例如，如果计算发现四声子贡献为负或物理上不合理，那模型很可能有问题。
    4. 与已知极限比较：计算出的热导率是否在理论上下限（如非晶极限、单声子气体模型预测）之内？温度依赖的指数是否合理？

6. 工作流总结与拓展应用

回顾整个研究，我们建立了一个用于强非谐性材料热输运研究的稳健工作流：第一性原理计算提供基准 → TDEP获取温度依赖的精确力常数 → 包含四声子散射的BTE计算热导率 → MLFF拓展至高压等极端条件 → 结合电子结构分析揭示微观机理。

这套方法的价值不仅在于解释了CuCl这一个特例，更在于其普适性。它适用于任何传统方法可能失效的强非谐性材料，例如：

• 热电材料：寻找本征低热导的材料，如SnSe、BiCuSeO等。
• 热障涂层材料：如稀土锆酸盐，其低热导往往与非谐性相关。
• 新型量子材料：某些拓扑材料或关联电子体系可能具有奇异的声子行为。

在实际操作中，计算资源的分配需要规划。TDEP+四声子散射的计算成本很高，通常需要在高性能计算集群上进行。一个实用的建议是：先小规模试算。用较小的超胞、较粗糙的q网格和较低的截断能做一个快速测试，评估四声子散射的贡献是否显著。如果贡献小于5%，或许可以忽略；如果像CuCl这样贡献超过80%，那就必须投入资源进行全精度计算。

最后，这项研究也展示了机器学习势函数在计算材料学中的强大潜力。它不再是黑箱，而是连接高精度量子计算与大规模有限温度/压力模拟的可靠桥梁。将MLFF与TDEP、BTE结合，为我们打开了一扇高效、精准研究复杂条件下材料热物性的大门。未来，结合主动学习策略，这套流程有望实现全自动的高通量筛选，加速下一代热管理材料和热电材料的发现。