机器学习势函数在暗物质探测中的应用：计算晶体缺陷存储能

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对接的量子机器学习势函数将是模拟极致低温下探测器响应的一个有前景的方向。这项工作展示了一个完整的闭环从一个具体的物理问题探测器能量响应出发驱动发展一种前沿的计算方法机器学习势函数通过大规模高精度模拟获得关键物理参数缺陷存储能分布最终反馈回实验用于修正物理分析并提出了新的检验实验。它体现了计算物理在现代基础科学实验中日益增长的作用——不再仅仅是解释实验而是开始预测和指导实验。