Study of Multiple Physical Constraints on the Duflo-Zuker Nuclear Mass Model

原子核质量(或称结合能)不仅是研究核内结构信息的有效探针，还是天体物理中研究重元素起源的关键物理量。一方面，核质量给出了核子关联、壳效应、形变等关键信息^[1]。另一方面，精确的核质量是核合成网格计算的重要输入^[2-4]，这对理解21世纪尚未解决的11个重大物理难题之三有着重要积极意义^[5]，即“宇宙中从铁到铀的元素是如何产生的”。因此，质量测量一直是核科学研究的热点。

实验上可通过直接或间接测量的方法获得约2 500个核素的质量，而不同理论预言给出在自然界中可存在的核素数目约是7 000~10 000^[4,6]。目前由于这些未知核素集中在质子滴线区、丰中子区和超重核区，目标核寿命短、产生截面小，精确测量其质量存在着许多技术上的困难。在可预见的将来，不稳定核素质量仍需依赖理论支持。

原子核质量的理论模型通常可分为两类。一类是全局型，另一类是局域型。全局型中精度较好的是宏观-微观和纯微观模型^[1]，宏观-微观模型的宏观项指的是描述结合能的平滑项，微观一般指向壳层和对能修正。纯微观模型的经典思想为给定核子相互作用势场，在非相对论下求解薛定谔方程获得结果，在相对论下求解狄拉克方程。全局型质量模型包括Duflo-Zuker(DZ)^[7]模型、有限力程小液滴模型^[8]和平均场理论^[9]等，其标准差约在300~600/keV^[10]。局域型则利用相邻核素间的质量关系间接获得未知核素质量，其理论精度约在150/keV以内，包括有著名的Garkev-Kelson(GK)关系^[11-12]、剩余中子-质子相互作用公式^[13-14]、同位旋多重态公式^[15-16]。国内课题组在局域质量关系中做出了许多重要的发现和拓展，包括对GK关系做出进一步推广并确定最优化外推路径的研究^[17-18]和镜像核外推公式的研究^[19-20]。关于局域质量关系的更多细节可查阅最新综述^[21]。虽然宏观-微观模型的物理思想广度和深度要明显高于局域型，但在一定质量数区域内局域型的结合能精度远远高于全局型。近年来，随着核物理和机器学习交叉学科的兴起，核质量模型的精度获得巨大改善^[22-31]，达到100/keV以内，即满足核合成计算所需的精度需求。但是，近几年一些科研人员提出了一个尖锐的问题：模型在已知区域核素的精度和未知区域的预言能力是否存在相关性？

不幸的是，Sobiczewsk等^[31-32]基于当前流行的数十个质量模型的系统性研究时发现模型在已知区域核素的精度和未知区域的预言能力的相关性较小^[33]。这在一定程度上表明单一追求结合能的精度可能并不会提升模型计算未知核的预言能力。另外，质量模型还存在着参数化模型的内部矛盾，即过拟合现象。比如，即使DZ-type质量模型计算结合能的精度较高，但质量相关量α衰变能的精度相对较低^[34]。不仅如此，不同质量模型的外推差异显著，在中子滴线附近差异能达到几十个MeV^[34-36]。这些现象在一定程度上表明了质量模型存在不可忽视的过拟合问题。虽然过拟合的模型能够完美重现拟合样本集的数据，但是在预测样本集外的数据时，误差将随着外推距离的增大而增大。

通常，缓和过拟合的主要方案是提供更多的拟合样本，如此可避免算法在某一区间内过度“学习”，进而达到自由参数全局最优而非局部最优。因此，我们在核质量领域中引入多目标优化算法，将α衰变能和GK关系作为两个重要条件去约束质量模型。本文旨在通过引入多个物理约束条件，提升质量模型在远距离外推时的预言能力。

Duflo-Zuker(DZ)系列是目前研究中使用率较高的质量模型，原因在于其具备较高的理论精度和外推性能，其计算成本也要远远地小于其他同等精度的质量模型。DZ系列模型包括DZ10、DZ28和DZ33等^[37-38]，本文采用DZ10模型为例，给出其多目标优化下的结果。

DZ模型由宏观部分的液滴项和微观壳模型 monopole哈密顿项组成。DZ10模型对一个原子核的基态结合能可以如下表示：

其中：$ a_{i}^{} $是由拟合确定的自由参数；$ V_{\rm{C}}^{} $、$ M + S $、${M}/{\rho }$、$ V_{\rm{P}}^{} $、$ V_{\rm{T}}^{} $、$ V_{{\rm{TS}}}^{} $分别代表液滴模型里的库仑能、体积能、表面能、对能、对称能和对称能修正能；$ s_3^{} $、${s_3^{}}/{\rho}$和$ s_4^{} $表示球形修正项；$ d_4^{} $是形变修正项。

人们在职场和生活中往往会面临许许多多的优化问题。比如，公司收入利润最大化、生活成本最小化等问题。然而，正如事物之间具有相关性，对某一目标要求收益最大化则必然会导致其他相关目标收益降低，因此需要对不同目标间的冲突关系做出阐述，这类同时考虑两个目标以上的优化问题称为多目标优化。多目标优化算法被广泛应用在经济学、工程学等学科^[39-41]。多目标优化算法定义了帕累托前沿解(Pareto front)，又被称为非劣解解集。非劣解解集中包含了一系列基于不同目标补偿关系所确定的解，非劣解的任何一个解几乎不可能使所有目标同时达到最优的效果，其核心思想是从各个冲突的目标函数条件下寻找尽可能多的补偿解，最终提供用户尽可能多的选择。不失一般性，一个多目标优化问题可以如下表示：

其中：$ f_{i}^{}(x) $是目标函数；$x = (x_{1}^{},x_{2}^{},\cdots ,x_{m}^{})$为决策变量；每个决策变量被限制在下界(lb)和上界(ub)的解空间($ \varOmega $)中。由于多目标优化中涉及多个目标函数的比较，因此无法像单目标优化中单独比较一个函数值就可确定个体的优劣。帕累托支配关系是多目标优化算法的核心要义，其定义解空间中个体的优劣关系。假设解空间中有两个个体$ s $和$ t $，个体$ s $支配(优于)个体$ t $, 并且仅当对$\forall i = 1,2,\cdots ,m$, $ f_i^{}(s)\leqslant f_i^{}(t) $，并且$\exists j \epsilon [ 1,2,\cdots ,m]$，使得$ f_j^{}(s) < f_j^{}(t) $。

目前，Non-dominated sorting genetic algorithm-II(NSGA-II)仍然是第二代进化算法中使用率最大、引用率最高的算法^[42-43]。因此，本文选择此算法对结合能和α衰变能进行多目标优化。α衰变能可由母核和子核的质量差给出，即$Q_{\alpha}^{} = BE(Z-2, N - 2)+ BE_{\alpha}^{} - BE(Z,N)$。目标函数的构造由理论值和实验值的标准差给出：

其中：$ \sigma_{K}^{} $为计算理论值和实验值的标准差公式；$ \hat{N} $为实验样本数，结合能数量为2 340个，偶偶核α衰变能数量为275个；$ K $是对结合能或其他物理量的代称。结合能和α衰变能的实验数据取自国内的原子核质量评估AME2020^[44]。图1为通过NSGA-II算法得出的DZ10模型非劣解结果图。横、纵轴分别代表结合能和α衰变能的标准差。首先，图1中的左上第一个解可理解成既是结合能精度最优解又是α衰变能精度最劣解，而右下的最后一个解则相反。我们可知质量模型中结合能和α衰变能的关系是相互冲突的，即结合能精度的提升是以劣化α衰变能精度为代价。这告诉我们以往的单目标算法在拟合结合能确定参数时，是以质量为最高优先级的局部最优化，而不是考虑结合能和α衰变能多个物理量的全局最优，这可能会使得过拟合问题加重，增大模型在外推时的不确定度。其次，在考虑了α衰变能的数据集后，可获得一系列计算α衰变能更占优的解。实际上，α衰变能是α衰变寿命的一个重要决定因素，以往的计算中发现α衰变能0.01/MeV的不确定度最终可影响α半衰期一个数量级，这些依赖模型本身获得的α衰变能占优解可以对预测重核和超重核的α衰变能提供一定支持。近年来，由于不同理论模型在预言实验值之外的区域时往往差异过大，因此理论学家则需要给出模型的系统和统计不确定度来对模型自身的不确定度给出定量描述。基于这几十个非劣解，我们可从给出模型的统计不确定度。表1为DZ10模型参数的不确定度。第二列为参数的统计不确定度；第三列为质量最优解的参数旨在给出对比。

基于非劣解解集的统计性，图2系统性地给出了新的核素边界、双中子分离能和非劣解的外推差异等物理信息。第一，液滴线可利用核子分离能小于零确定，由于偶偶核往往比奇A核稳定，所以通过双核子分离能确定的液滴线包含了通过单核子分离能确定的液滴线。图中，$ P(S_{2{\rm{n}}}^{}<0) $为通过双中子分离能确定的该核素是否为非束缚核的概率，$ P(S_{2{\rm{P}}}^{}<0) $为通过双质子分离能确定的该核素是否为非束缚核的概率。第二，我们系统地计算了非劣解的双中子分离能，并取加权平均，大小通过图中左侧数值颜色棒表示。第三，我们采用非劣解解集中的质量最优解和次质量最优解对应参数组，即图1中最上方相邻的两个解，并计算其之间的的结合能差异，在右下子图表示。此子图清晰表明，在有实验值和丰质子区时，外推的差异不超过$ 0.5 $/MeV。但是，在超重核区和丰中子区时，差异则会明显增大。我们考虑到这两个解具备明显偏向质量精度的特征，因此可能会导致大规模外推时过拟合现象加重，最终增大模型外推的差异。于是，我们从图1的中间区域任取两个解计算外推的差异。通过比对图2中子图和图3，可知中间区域的外推差异要明显地小于质量优解组的外推差异，表明利用平衡结合能和α衰变能的非劣解时外推时的不确定度较小。

在前文讨论的基础上, 可以对质量模型的约束条件做进一步的探索, 即通过引入局域质量公式来约束模型。本文将Garvey-Kelson关系作为限制核质量模型的另外一个条件。过去，Barea等^[45]就提出可对全局型模型拟合时附以局域型约束条件，以此增强全局型模型计算相邻核的连贯性和外推时的预测性，但文中讨论的实践方法较粗糙且冗余。之后，全局型模型往往将GK关系作为后验样本而非限制条件。我们借助多目标优化算法，进一步将Garvey-Kelson关系作为一个严格的限制条件去优化宏观质量模型。本文具体采用的GK关系如下所示：

其中：$ N $和$ Z $分别是中子数和质子数；$ M $是质量；$ \varDelta M_{{\rm{GK}}}^{} $是相邻6个核素质量加加减减抵消后的不确定常数，通常在$ \pm $150/keV以内。

$ \varDelta M_{{\rm{GK}}}^{{\rm{expt}}.} $是从AME2020中将6个实验已知的核素质量带入到式(4)中所得到的GK实验观测量，$ \varDelta M_{{\rm{GK}}}^{{\rm{theo}}.} $是DZ10模型计算与实验值相对应的GK理论观测量。我们将此约束条件式(5)作为通过GK关系限制DZ质量模型的第三个目标。实际中，基于剩余核子相互作用相互抵消假设，衍生出了许多版本的GK关系，可多达数十个。本文简单采用了式(4)将其作为可实现案例，表明了利用GK关系多目标限制质量模型的可行性，当然也可将式(4)替换成其他版本的GK关系，这还有待未来做出进一步的检验。图4为多目标优化算法下以结合能、α衰变能和GK关系给出的非劣解解集图。青色方框包裹的为质量最优解。上文已经谈论过多个物理量的平衡解在外推上可能更可靠，但此说明还仍待进一步确认，下文将对此关键问题做出细致分析。面对如此多的非劣解，如何从中寻找出更有物理意义、更能给出可靠外推值的解是十分关键的问题。由于丰中子核素具备较大的同位旋不对称度，对理论和实验而言预测和测量这些核素质量都是巨大的挑战，且不同理论模型的差异往往集中于丰中子区域。

首先，我们从近几年的实验中提取出尚未包括在AME2020中的丰中子核质量数据^[46-48]。然后采用DZ10的质量最优解计算这些核的标准差并将其作为一个理论基准值，若其余的非劣解给出比此基准值精度更好的标准差时，将保存这些非劣解并归入到潜力解解集。事实上，非劣解的数目高达几百个，但筛选出来的只有几十个潜力解，这些潜力解自然可以提供一些有用信息。表2从局部和整体两方面对潜力解给出系统性分析结果，第二至第四行针对局部丰中子核数据，最后三行针对整体拟合样本数据。其中：第二行是质量最优解针对丰中子核计算得出的标准差；第三行是潜力解解集中所能给出的最好标准差；第四行是所有潜力解结合能标准差的统计分析，其中第一项(0.304/MeV)为标准差的平均值，第二项(0.095/MeV)是标准差描述其浮动程度；第五、六、和七行分别是潜力解解集计算拟合样本中结合能、α衰变能和GK关系的统计分析。首先，最好潜力解计算的丰中子核质量标准差远远高于质量最优解结果，约提升70%的精度。其次，观察表中第四行可知，这些数据点的变化浮动较小，且潜力解解集对α衰变能和GK关系目标的平均值集中于各自的中间区域，即平衡区域。这在一定程度上验证了上文的论点，平衡多个物理量下的全局解更有可能给出稳定可靠的外推值。

上文在考虑结合能和α衰变能的关系下，采用了平衡两者关系的两个全局最优解计算其之间的结合能差异，发现结合能在外推时的差异已经减小，在某种程度上表明过拟合得到一定改善。自然，我们需要考虑在加入了结合能、α衰变能和GK关系下，即更多物理数据的约束下，提取的非劣解的结合能差异情况。图3(b)为从潜力解解集(即中间区域)任挑两组解计算的结合能差异图。相比于图3(a), 图3(b)在外推时的差异更小，表明加入GK下的可行性，降低了外推的不确定度。

超重核(SHN)的合成一直是核物理的前沿热点研究之一。研究超重核的途径主要是通过α衰变，原因在于α衰变是超重核的主要衰变模式。我们基于潜力解解集对超重核($ Z \geqslant 104 $)的α衰变能进行了理论计算。图5为潜力解解集对所有超重核计算α衰变能标准差分布图。由分布情况可知，此图粗略符合高斯分布，质量最优解处于分布中后端的位置，而大部分的潜力解集中在质量最优解的左侧，少部分的潜力解在质量最优解的右侧。这表明在大概率下潜力解所计算超重核的α衰变能会比质量最优解更加精确。另外值得一提的是，原子核的壳演化一直是国际上的重点研究方向，研究表明远离$ \beta $稳定线的弱束缚核素可能会出现壳层的演化。DZ质量模型涉及到关于壳层结构的复杂价核子计算，提前将Z(N) = 8, 20, 28, 50, 82, $ N $= 126, 184的壳层信息输入到质量模型的价核子计算之中，而实验中出现的新幻数尚未被考虑到DZ模型的计算中。本文也对此研 究进行了初步的尝试。我们在DZ10的壳层中添加了$ N = 152 $的中子壳，并对$ 126\leqslant N \leqslant 184 $的α衰变能进行计算，所得结果和实验值差别相对大些，此方法还有待未来进一步检验。

科学家们对铁以下元素合成的机制和天体环境已经有了基本共识，然而对于比铁重的元素在宇宙中是如何合成的，还存在着许多问题。快中子俘获过程(r-process)被认为是产生超铁元素的重要途经之一。发生r-process的天体环境被认为是在中子星的合并或是超新星的爆发之中，但是这两种天体环境都存在一定的局限性。另外，r-process的模拟需要大量核物理数据，且所需核素集中在极端丰中子区，然而实验上这些核的基本物理量比如质量、反应率数据极度缺乏，所以必须要依赖于理论模型的预测。

本文选用的r-process场所是超新星爆发产生的中微子驱动星风^[49]，使用的核合成模拟计算代码是NucNet^[50]。其中初始温度设置为$ T = 10 \; {\rm{GK}} $，密度随时间变化的表达式为

其中：$\rho_0^{} = 1.498\,5\times 10^6$，$ \rho_1^{} = 1.5\times 10^3 $，$ \tau_0^{} = 0.035 $，$ \tau_1^{} = 1 $为给定初始参数，单位为g/cm³。

在核物理输入方面，我们采用的核反应截面数据均来自JINA REACLIB数据库^[51]。然而对于原子核的质量，我们并没有采用JINA REACLIB里的实验值加有限力程小液滴模型理论值结合的质量表，对实验尚未测得的核，我们将采用上文得到的DZ10非劣解解集去计算理论质量值，然后运行r-process模拟。根据比较丰度理论值和实验值之间的吻合程度，我们则可进一步地挑选出在外推上更可靠的理论质量表。比较图6中丰度的理论值和实验值可知，DZ10的非劣解已经基本重现了丰度实验值的分布，尤其是再现了$ A = 130 $和$ A = 195 $的两个峰。通过比较理论值和实验值的吻合情况，我们可以进一步地挑选出质量预言更可靠的非劣解。

本文旨在通过添加严格的物理约束，减轻核质量模型的过拟合程度，从而降低模型远距离外推时巨大的不确定度。在多目标优化算法的框架下，引入α衰变能和局域质量公式Garvey-Kelson关系，起到大幅提升α衰变能的精度，增强模型外推能力的效果。此外新核素边界和模型统计不确定度通过多目标方法的统计特征给出。我们从海量的非劣解解集中筛选出更具外推能力的解(称为潜力解)，并对潜力解做出严格的统计分析，我们得出平衡多个物理量的全局解比单一的质量最优解的预言能力更高。这些物理限制的引入在一定程度上减轻了质量模型的过拟合问题，因此提升了模型在外推时的可靠性。这将会对丰中子核的可靠预言提供一定理论支持。另外，我们还通过r-process的模拟计算更进一步地挑选出质量预言更可靠的非劣解。