一种从热容数据中提取玻璃振动态密度的稳健方法：在非晶态二氧化硅中的应用

《Computer Physics Communications》：A Robust Method for Extracting the Vibrational Density of States of Glasses from Heat Capacity Data: Application to Vitreous Silica

　　

• •
  从 Cp 数据中提取 VDOS 的自动正则化参数选择。
• •
  Bootstrap 重抽样（B=200）得出首个严格的 VDOS 置信区间。
• •
  系统比较一阶/二阶 Tikhonov 正则化与 Surovtsev 方法。
• •
  与非弹性中子散射数据验证一致，相关系数 r ≈ 0.98。
• •
  提供开源 Python 脚本以实现完全可重复性。

引言

振动态密度（VDOS）是表征凝聚态系统动力学性质的最基本物理量之一[[1], [2], [3]]。对于非晶材料，缺乏长程结构有序性使得标准晶格动力学方法不再适用，需要通过非弹性中子散射或 X 射线散射 [4,5] 进行实验表征，或从热力学测量中提取。玻璃在低温下表现出独特的异常现象。在毫开尔文温度区间，异常热性质通常用隧道态来描述 [6,7]。在较高的低温范围内，低频振动模态相对于 Debye 预言的 g(E) ∝ E² [8] 存在过剩，由此产生了玻色子峰区域。这种过剩表现为约化态密度 g(E)/E² 在数 meV 能量处出现一个宽峰 [9,10]。对其位置、强度和谱形的定量表征为区分不同竞争物理描述提供了有用的约束条件[[11], [12], [13]]。

一种有吸引力的互补方法利用 VDOS 与热容之间的热力学关系。在谐运动近似下，定容热容 C_v(T) 是 g(E) 的积分变换，而本文所考虑的低温范围内 C_p 被用作 C_v 的近似。精密量热法比专用散射设施具有更大的可获得性，但从 C_p(T) 反演 g(E) 是一类 Fredholm 积分方程，因此是不适定问题 [14,15]。数据微小扰动可能导致无约束谱分量出现巨大变化，因此需要正则化。Tikhonov 正则化提供了一个自然的框架 [16,17]，Surovtsev 曾为非晶二氧化硅建立了一种有影响力的一阶实现方案，其中正则化强度通过联合考察重建诊断和振荡诊断来确定 [18]。该文献中引用的约 3% 精度指的是其均方根相对重建度量 χ_s，这与下文使用的平均绝对相对误差是不同的指标。

因此，核心计算问题不仅仅是最小化前向残差，而是要识别在物理合理的正则化选择和数值实现下保持稳定的谱观测量。本工作通过以下四个互补要素来解决这一问题：(i) 直接在增广系统上求解的非负 Tikhonov 正则化公式；(ii) 一阶和二阶正则化与 Surovtsev 型参考分支的比较；(iii) L 曲线、预测交叉验证和固定参考重抽样诊断 [19,20]，以区分参数选择敏感性与重抽样稳定性；(iv) 通过与独立获得的绝对 VDOS 测量进行定量验证，该测量采用非弹性 X 射线散射核共振分析（NRAIXS）方法获得 [21]。我们提供了一个带有完整文档的 Python 软件包，以复现本文报告的主要量热提取结果、参考输出、L 曲线分析以及交错预测交叉验证诊断。

章节摘要

前向问题：从 VDOS 到热容

出发点是在谐运动近似下振动态密度 g(E) 与低温热容之间已确立的关系 [8,22,23]。对于量子谐振子系综，配分函数给出了能量为 E 的模式的贡献。在本文使用的温度范围内，C_p 被用作 C_v 的近似 [24]，这与整个分析中采用的量热公式一致。基本关系式为：

热容数据

非晶二氧化硅的热容数据来源于 Flubacher 等人的经典测量 [25]。提供的数据集包含 33 个测量点，温度范围为 2.344 至 18.98 K。以 cal/(mol·K) 为单位的表格热容值通过 4.184 J/cal 进行转换，再除以三，以得到三原子 SiO₂ 化学式单元的单原子值。所有参考计算均使用相同的 33 个数据点。

用于独立验证，我们采用 Chumakov 等人报道的环境条件非晶二氧化硅约化 VDOS。

正则化诊断与参考强度

图 2a 展示了一阶和二阶正则化的约束平均相对重建误差 ε(α)。在考察的 α 低值范围内，约束曲线未表现出明确定义的内部极小值。相反，ε 接近约 0.256% 的近似平坦底值，仅在正则化足够强以致恶化量热重建时才会上升。因此，仅靠前向拟合误差在较宽的 α 区间内几乎不提供区分能力。

与 Surovtsev 方法的关系

本文使用的 Surovtsev 型分支不应被解释为对文献 [18] 的精确复现。Surovtsev 采用了一阶光滑泛函，在没有显式非负性约束的情况下求解了所得线性方程组，并通过目视平衡 χ_s 和 γ 选择了 α ≈ 3 × 10^?5。在本文的比较中，使用相同的 α 数值在共用的约束增广系统求解器上进行评估，使得三个参考分支之间的差异反映正则化方案本身的差异而非求解器的差异。

结论

我们开发了一个受约束且可复现的框架，用于从低温热容测量中提取低能振动态密度。主要反演过程在增广非负最小二乘求解中强制施加非负性约束，使用原始有限差分 Tikhonov 惩罚项，并将科学点估计与固定 α 的重抽样统计量分离。

对于非晶二氧化硅，一阶、二阶以及 Surovtsev 型参考分支均能复现……

补充信息

附录 A 记录了数值公式、正则化诊断、条件数分析、数值路径敏感性以及独立的 NRAIXS 验证。补充材料包括一个基于热容的 VDOS 提取通用 Python 软件包。该软件包包含两个脚本，仅产生数值 CSV 输出；不生成图形。

vdos_extraction.py：主要的一阶和二阶约束 Tikhonov 反演。该脚本接受 CSV 或 Excel 格式的热容……

CRediT 作者贡献声明

Omar Benzine：撰写——初稿，可视化，软件，方法论，调查研究，形式化分析，数据管理，概念构思。