机器学习势函数在暗物质探测中的应用：计算晶体缺陷存储能
1. 项目概述：当机器学习势函数遇上暗物质探测
在粒子物理与凝聚态物理的交叉前沿，有一个看似微小却至关重要的物理细节，正困扰着新一代的暗物质与中微子探测实验：当一个来自宇宙的弱相互作用粒子（WIMP）或一个反应堆中微子，与探测器晶体中的一个原子核发生碰撞时，这个原子核获得的动能（通常在10-100电子伏特量级）究竟有多少会转化为可被测量的热信号（声子）？又有多少能量会“悄无声息”地存储在晶体内部新产生的缺陷中，从而“丢失”了？这个“能量存储”效应，直接关系到探测器能量标定的准确性，进而影响对暗物质粒子质量、相互作用截面等关键物理参数的提取。
传统上，回答这个问题依赖于分子动力学模拟。但这里有个两难困境：使用经验势函数（如Tersoff、EAM）速度快，但难以准确描述复杂氧化物（如钨酸钙，CaWO4）中多体相互作用的细节；使用第一性原理密度泛函理论计算精度高，但其巨大的计算成本使得模拟包含数万个原子、持续数皮秒的位移级联过程成为不可能的任务。这正是我们这项工作的切入点：我们开发并应用了一种基于线性机器学习模型的原子间势函数，对CaWO4晶体中由100 eV量级核反冲引发的位移级联进行了大规模、高精度的分子动力学模拟。我们的目标很明确——定量计算晶体缺陷存储的能量，并评估其对暗物质和中微子反冲能谱的扭曲效应，为下一代高精度实验提供不可或缺的理论修正。
2. 核心思路与技术路线拆解
2.1 问题定义与物理图像
让我们先构建一个清晰的物理图像。在低温CaWO4晶体探测器中，一个入射的暗物质粒子或中微子与钨（W）原子核发生弹性散射。钨原子核获得动能，成为初级撞出原子。这个被撞出的原子在晶格中穿行，与沿途的其他原子发生一系列碰撞，形成一个“位移级联”。这个过程就像在平静的池塘里投入一颗石子，激起的涟漪会扰动周围的水分子。在晶体中，这种扰动可能导致原子离开其平衡位置（产生空位），或挤入晶格间隙（产生间隙原子），形成所谓的弗伦克尔缺陷对。
关键在于，创造这些缺陷需要消耗能量。这部分能量来自初始的核反冲动能，但它并没有转化为晶格振动（热能/声子），而是以缺陷的势能形式“存储”在了晶体中。因此，探测器实际测量到的“可见能量”等于初始反冲能量减去存储于缺陷中的能量。如果忽略这一修正，实验测得的能谱峰值将发生偏移，能谱形状也会发生扭曲，导致对物理参数的错误解读。
2.2 技术路线：从第一性原理到大规模分子动力学
我们的技术路线是一个典型的“自上而下”的多尺度模拟流程，其核心是用机器学习势函数这座“桥梁”，连接高精度但小尺度的第一性原理计算与大规模但需快速势函数的分子动力学模拟。
第一步：构建第一性原理数据库一切高精度模拟的基石是可靠的数据。我们使用VASP软件包，在广义梯度近似（GGA-PBE）泛函下，对CaWO4晶体进行了系统的密度泛函理论计算。为什么选择PBE？虽然更精确的杂化泛函（如HSE06）存在，但其计算成本高出1-2个数量级。对于需要构建包含数百个构型的训练数据库的任务，PBE在精度与效率之间提供了最佳平衡。我们后续也用HSE06检查了关键缺陷的形成能，确认了PBE结果的可靠性，这步交叉验证至关重要。
数据库的构建并非随机采样，而是有针对性地包含了与位移级联相关的关键原子构型：
- 弗伦克尔缺陷对：这是低能反冲最可能产生的主要缺陷类型。我们系统搜索了钨、钙、氧原子的间隙位置，计算了其与相应空位配对的形成能。表I中的数据是后续评估机器学习势函数精度的黄金标准。
- 高温分子动力学快照：从完美晶体和含一个钨弗伦克尔对晶体的2000 K MD模拟中抽取构型。这确保了势函数能描述有限温度下的原子振动和动态行为。
- 近距离原子对构型：在模拟高能碰撞的初始阶段，原子可能被挤压到极近的距离。我们通过固定原子间距的约束弛豫，生成了这类极端构型，为势函数描述短程强排斥作用提供数据。
最终，我们得到了一个包含348个构型的数据库，每个构型都有对应的总能、每个原子上的受力以及晶胞的应力张量。这三类信息共同作为训练目标，确保学习到的势函数不仅能预测能量，还能准确预测原子受力和应力，这对于动力学模拟的稳定性至关重要。
第二步：设计线性机器学习势函数我们采用了线性机器学习模型，其核心思想是将每个原子的能量表达为其局部原子环境描述符的线性组合。模型的优势在于其简洁、高效且不易过拟合。这里有两个技术关键点：
混合描述符：我们使用了两种描述符的直和。
- 双体径向描述符（za）：快速计算，主要描述原子对之间的径向距离信息。
- 多体描述符（xa）：基于SO(4)双谱系，这是一种对原子局部环境旋转、平移和置换对称性具有不变性的数学描述。它将原子邻居的化学种类和空间分布信息编码进一个固定长度的向量中。我们设置了4.8 Å的截断半径，并限制角动量最大值为4，在表达能力和计算成本间取得平衡。
短程排斥势的平滑衔接：在位移级联的初始瞬间，原子间距可能小于1 Å，此时电子云重叠导致的排斥力非常强。标准的DFT计算（如VASP）对此描述并不完善。因此，我们引入了一个关键技巧：当原子间距小于1.3 Å时，逐渐关闭多体描述符的贡献；当间距小于1 Å时，用经典的Ziegler-Biersack-Littmark势来描述这种核-核之间的屏蔽库仑排斥。在1 Å到1.3 Å之间，用一个缓冲函数实现ZBL势与我们的2体核势之间的平滑过渡，确保势能面及其一阶（力）、二阶（声子）导数连续。没有这个处理，模拟高能碰撞事件时原子可能会非物理地“穿过”彼此。
使用MILADY软件包进行线性回归拟合后，我们的势函数在测试集上达到了平均能量误差0.474 eV，平均原子力误差0.131 eV/Å的精度。对于涉及数百个原子、总能量变化在数千eV的模拟来说，这个精度足以可靠地分辨出10 eV量级的缺陷存储能差异。
第三步：大规模分子动力学模拟与采样有了可靠的势函数，我们使用LAMMPS分子动力学软件（通过MILADY插件调用我们的势函数）进行位移级联模拟。
- 模拟体系：构建了一个包含49152个原子（约8.5×8.5×9.2 nm³）的近立方晶胞。这个尺寸确保了位移级联产生的热峰和应力场不会影响到周期性镜像，同时缺陷有足够的空间弛豫。
- 温度设定：探测器工作在约10 mK的极低温下，但原子并非静止，而是具有零点振动能。为了在经典的分子动力学框架中近似考虑量子效应，我们根据一个经验规则，将等效经典温度设为CaWO4德拜温度（354 K）的3/8，即133 K。在这个温度下初始化体系并进行热化，使得原子的均方位移与量子基态近似匹配。
- 反冲模拟：从热化平衡的体系中，选取中心的一个钨原子作为初级撞出原子，赋予其特定的动能（20, 40, 60, 81, 112, 160 eV）和方向。然后让体系自由演化8皮秒，使用变时间步长算法，并以最大原子位移0.01 Å作为收敛判据。模拟结束后，将体系“淬火”到局部能量极小点，通过计算与完美晶体基态的能量差，得到存储在缺陷中的能量。
- 方向与初始构型采样：钨原子位点具有4度旋转对称性，因此我们只需在1/4的立体角内均匀采样32个反冲方向即可覆盖所有可能性（见图1左）。此外，为了平均掉热涨落的影响，对每个反冲方向和能量，我们从热化轨迹的最后10 ps中抽取了多个（通常为10个）不同的初始原子速度构型进行独立模拟。这种充分的采样是获得可靠统计分布的前提。
3. 模拟结果与物理发现
3.1 缺陷存储能的统计分布
通过对大量模拟轨迹的分析，我们得到了不同初始反冲能量下，缺陷存储能的概率分布（图2）。几个关键发现颠覆了此前一些简化假设：
平均存储比例约为10%：对于我们所关心的能量区间（40-160 eV），平均约有10%的初始反冲能量被存储在晶体缺陷中。这意味着一个标称100 eV的反冲事件，探测器实际探测到的“可见能量”只有约90 eV。这是一个必须进行校正的系统性偏移。
存储能量与反冲能量呈非线性关系：如果存储机制是简单的比例关系，那么160 eV反冲的存储能应该是81 eV的两倍左右。但我们的模拟显示，81 eV反冲平均存储7.5 eV，而160 eV反冲平均存储11.0 eV，后者并非前者的两倍。这揭示了“位移阈值能量”的存在：只有当反冲原子传递给晶格原子的能量超过某个阈值（足以将其踢出晶格位置）时，才会产生缺陷。随着反冲能量增加，产生缺陷的效率（单位动能产生缺陷的几率）并非恒定，导致了这种非线性。
分布平滑，非离散化：存储能的分布是连续且平滑的，而非集中在几个离散的值上。这否定了“只有少数几种特定缺陷构型”的简单图像。通过对模拟末态原子构型的直接观察，我们发现反冲轨迹周围区域出现了局部的非晶化，即原子排列长程无序。这表明低能位移级联产生的不是几个孤立的点缺陷，而是一个小的、紊乱的原子团簇，其精确构型和能量因碰撞细节（方向、初始热振动相位）而异，从而产生了连续的存储能分布。
方向依赖性不明显：图1右图展示了112 eV反冲在所有方向上的存储能分布（经立体角加权平均后）。令人惊讶的是，存储能与反冲方向没有表现出强烈的相关性。尽管钨原子在CaWO4晶格中的局部环境是各向异性的，但在我们模拟的能量尺度下，反冲过程似乎足够剧烈，以至于方向性的细节被碰撞的随机性所掩盖。这简化了后续的能谱修正，因为我们不需要引入复杂的角度依赖修正因子。
3.2 从离散模拟到连续能谱修正
我们的MD模拟只针对有限的几个离散反冲能量（20, 40, 60, 81, 112, 160 eV）进行。然而，暗物质和中微子的反冲能谱是连续的。为了将我们的结果应用于整个能谱，我们需要一个插值方案。
我们的方法是：
- 对每个模拟能量Ei，用参数化模型拟合存储能分布。我们尝试了对称高斯、非对称高斯等多种模型，发现模型的选择对最终卷积结果影响甚微，这说明了我们基础数据的稳健性。
- 将模型参数表达为初始反冲能量E的连续函数。我们尝试了多项式、高斯函数等多种拟合形式。
- 引入缺陷产生概率pd(E)。我们发现，pd(E)从低能下的0平滑地过渡到高能下的1，其形状可以用一个费米函数完美拟合（图3）。这个函数本质上描述了“超过位移阈值的概率”随能量升高而增加的过程。
最终，对于任意初始反冲能量E，我们都有一个对应的、归一化的存储能概率分布函数P(S|E)，其中S是存储能。实验观测到的“可见能量”E_vis的谱形，需要对原始理论反冲能谱dR/dE进行如下卷积计算： dR/dE_vis = ∫ dE (dR/dE) * ∫ dS P(S|E) * δ(E_vis - (E - S)) 这相当于用我们的存储能分布“模糊化”了原始的理想能谱。
4. 对物理实验的直接影响与验证途径
4.1 对暗物质与中微子探测能谱的修正
我们将上述修正流程应用于两个典型的物理过程：
2 GeV质量暗物质粒子的弹性散射谱：假设标准晕模型和自旋无关相互作用。修正后的能谱（图5黑色实线）相对于未修正的谱（黑色虚线）发生了明显的形状扭曲。在20 eV阈值以上，修正因子在0.9到1.1之间变化，即产生了高达±10%的扭曲。这种扭曲与某些新物理场景（如中微子非零磁矩效应）可能产生的信号畸变在形态上具有相似性。如果不进行本底扣除或联合分析，这种源于探测器材料的效应可能会被误认为是新物理的信号，或掩盖真实的新物理信号。
反应堆中微子相干弹性散射谱：计算结果（图5红色曲线）显示，其修正幅度和形状与2 GeV暗物质的情况非常相似。这对于像NUCLEUS这样旨在以高精度测量CEνNS的实验至关重要。我们的计算为从实验数据中提取精确的弱混合角或寻找超越标准模型效应，提供了一个必须考虑的实验室本底修正。
4.2 一个独特的实验验证方案：中子俘获核反冲
理论计算需要实验验证。我们指出，利用热中子俘获过程产生的单能核反冲，可以为我们的模拟提供理想的检验平台。具体原理是：钨-183原子核捕获一个热中子后，形成激发态的钨-184，随后退激放出多个γ光子。根据动量守恒，发射γ光子会导致钨原子核本身受到反冲。通过测量不同γ射线级联的角关联，可以精确推算出最终钨核的反冲能量，产生81、112和160 eV等离散的单能反冲峰。
我们模拟了在一个0.75克CaWO4探测器、能量分辨率为5 eV的理想情况下，中子俘获反冲能谱的修正效果（图4）。可以清晰地看到，三个校准峰都向低能方向发生了大约10%的偏移。更重要的是，81 eV和160 eV峰的偏移量并非严格的2倍关系，这正体现了我们预测的非线性效应。
目前，基于此原理的CRAB校准方法正在推进实验。计划在维也纳TRIGA反应堆上，将低温稀释制冷机与热中子束流结合，进行测量。预计在几天的数据采集后，峰位的统计精度可以达到亚电子伏特级别。这不仅能以前所未有的精度验证我们的缺陷存储能计算，还能直接探测探测器在极低能区（~100 eV）的响应非线性。如果实验能进一步通过符合测量反冲钨核与特定γ射线的方向，甚至可能研究存储能是否具有微弱的各向异性，这将是检验我们模拟细节的终极测试。
5. 实操启示与未来展望
5.1 机器学习势函数开发的经验与教训
- 数据库的质量决定势函数的上限：构建训练数据库时，必须紧密围绕目标应用场景。我们的数据库不仅包含平衡态性质（晶格常数、弹性常数），更重点包含了非平衡态、高能过程相关的构型（弗伦克尔对、近距离原子对、高温快照）。如果只使用平衡态数据训练，得到的势函数几乎肯定无法正确模拟位移级联。
- 描述符与回归模型的选择需要权衡：线性模型配合精心设计的描述符（如双谱系），在保证精度的同时，计算效率远高于复杂的深度神经网络模型。对于需要运行成千上万次MD模拟的统计采样任务，计算速度是关键。我们的线性势函数在保持0.1 eV/Å力精度的同时，其计算速度比在线DFT计算快5个数量级以上，使得大规模采样成为可能。
- 短程相互作用处理是高压/高能模拟的生死线：对于涉及原子近距离碰撞的模拟，必须妥善处理短程强排斥势。简单地外推由中程数据训练出的势函数是危险的。引入ZBL势并实现平滑衔接，是一个经过验证的可靠方案。
- 验证、验证、再验证：势函数训练完成后，必须在一系列DFT未训练过的性质上进行测试。我们对比了PBE和HSE06计算的缺陷形成能，验证了势函数在关键性质上的可靠性。此外，还应测试声子谱、热膨胀系数等，确保势函数在应用温度范围内是可靠的。
5.2 对极低能粒子探测实验的启示
- “能量标定”需要重新审视：传统的探测器能量标定通常使用单能X射线或γ射线，它们通过光电效应沉积能量，主要产生电子空穴对或声子，不涉及核反冲和晶格位移。而暗物质和中微子信号是核反冲。我们的工作明确指出，核反冲的“可见能量”与沉积能量之间存在一个与材料、能量相关的“淬灭因子”，这个因子不仅包含传统的电子/声子转换效率，还应包含缺陷存储能这一项。未来实验的蒙特卡罗模拟和数据分析管线必须集成这一修正。
- 材料筛选的新维度：在设计和筛选下一代探测器的靶材时，除了高原子序数（提高相干散射截面）、低本底、易生长等传统指标外，低缺陷存储能或可预测的缺陷存储能也应成为一个考量因素。通过高通量的机器学习势函数计算，可以快速筛查不同晶体材料（如Si, Ge, Al2O3, CaWO4等）的缺陷存储特性，为实验选材提供理论指导。
- 面向更高精度的模拟需求：我们的工作基于经典的分子动力学，并通过设置等效温度来近似量子零点振动。对于更低温（~10 mK）的探测器，量子效应可能更加显著。路径积分分子动力学可以严格处理量子核效应，但其计算成本极高。发展能与PIMD对接的量子机器学习势函数，将是模拟极致低温下探测器响应的一个有前景的方向。
这项工作展示了一个完整的闭环：从一个具体的物理问题（探测器能量响应）出发，驱动发展一种前沿的计算方法（机器学习势函数），通过大规模高精度模拟获得关键物理参数（缺陷存储能分布），最终反馈回实验，用于修正物理分析并提出了新的检验实验。它体现了计算物理在现代基础科学实验中日益增长的作用——不再仅仅是解释实验，而是开始预测和指导实验。