基于协变能量密度泛函的250Cf 诱发裂变动力学研究

李泽宇

(西南大学物理科学与技术学院，重庆 400715)

对重核裂变的微观描述是低能理论核物理中最复杂的问题之一[1，2]。由于重核中通常包含着数以百计的核子通过多体相互作用耦合在一起，使用从头算方法(Ab initio)或者相互作用壳模型(Interacting Shell-Model)计算如此复杂的多体问题都显得捉襟见肘。因此，现代微观方法对于重核的研究通常基于原子核能量密度泛函(NEDFs)的框架。其中原子核密度泛函理论(DFT)和时间依赖的密度泛函理论(TDDFT)使人们能够完成对重核的静态性质以及动力学演化的自洽计算[3-11]，但经典的含时密度泛函理论计算只能给出原子核的最可几演化路径，并不能给出实验观测到的裂变产物呈现类高斯分布的结果。

一种可靠的恢复裂变动力学量子性的方法是时间依赖的生成坐标方法(TDGCM)。利用这种方法，将原子核波函数视为集体形变空间中内禀态的叠加，通过对叠加系数的动态演化即可再现裂变过程中的量子效应。基于高斯重叠近似(Gaussian overlap approxima)，从含时的Hill-Wheeler 方程出发可以推导出集体空间中局域的含时类薛定谔方程。求解此方程即可得展开系数的时间依赖关系，此时，裂变系统的动力学演化本质上取决于集体坐标，能量密度泛函，对相互作用的选择以及集体质量的计算[12]。

在本文中，我们利用协变密度泛函理论计算了250Cf 原子核的静态位能曲面，利用TDGCM+GOA 研究了集体空间中波函数的演化。在第二节概述了于计算裂变系统位能曲面、集体惯性和含时演化的理论模型。在第三节中，展示并探讨了对于250Cf 裂变的动力学计算结果。其中裂变碎片的质量分布呈现四峰结构。第四节中对研究结果进行了总结。

1 理论框架

核裂变可以被近似为一个只由几个集体自由度决定的缓慢绝热过程。通过绝热近似，原子核的内禀核子自由度与集体形变自由度解耦，这意味着我们可以将计算分为独立的两步。第一步，通过施加约束的方法计算原子核在特定形变下的内禀组态，以得到原子核的位能曲面以及剪裂线，子核核子数等集体参数。第二步，通过求解类薛定谔方程，计算集体波函数在集体空间中的动力学演化。在本文中，第一步通过基于PC-PK1[13]的协变密度泛函理论(CDFT)程序计算完成。第二步则通过Felix2.0 程序[14]计算模拟。

1.1 集体参数的求解

使用基于PC-PK1 的协变密度泛函理论，求解单核子Dirac 方程并结合BCS 方程便可得到单核子波函数，占据几率，体系能量等静态参数[13]。本文将Dirac 旋量φk(r)在各项异性的形变谐振子基上展开，通过自洽迭代求解。为了获得原子核的裂变位能曲面，需要在计算中对四极矩和八极矩施加约束。

其中＜H＞是总能量，＜Qk＞表示质量四极和八极算符的期望值，qk是对应的约束值，Ck是约束强度。实际计算中，定义形变参数β2和β3替代原本的四极与八极算符期望值：

其中R0=r0A1/3=1.2*A1/3(fm)。

为了实现类薛定谔方程的计算，需要计算原子核微观集体质量张量Bkl以及集体势场Vcoll。通过微扰推转近似(perturbative cranking approximation)，质量张量可以表示为[15，16]：

其中：

求和遍历所有的核子准粒子态。集体势场则通过CDFT计算结果减去零点能获得[17，18]。

1.2 含时类薛定谔方程

以原子核形变参数β2，β3为生成坐标，采用高斯重叠近似，从含时的Hill-Wheeler 方程出发可以推导出集体空间中局域的含时类薛定谔方程。

其中g(β2，β3，t)表示在(β2，β3)集体形变空间中的集体波函数。此方程的输入参数皆由相对论能量密度泛函自洽计算而得，求解此方程即可得集体波函数随时间演化的关系。由于位能曲面计算采用变分原理，整个曲面将被分为两部分。当目标约束形变相对较小时，计算结果为正常裂变所经过的组态。而当约束计算极端形变原子核时，将会得到原子核聚变道的结果[19]。这将导致最终得到的位能曲面上存在断崖式的能量骤降区域。

实际计算表明，由裂变道到聚变道的突变通常在脖子数为3 核子左右时发生。本文中，定义脖子数为3 核子的组态为断点，每一个断点都对应着一种断裂组态。断点所连成的线为剪裂线，对于剪裂线上任一面元ε，累计计算波函数的通量：

即可得该面元对应断裂组态发生的概率，考虑整个剪裂线即得最终产物分布。集体波函数的概率流可以通过连续性关系求得：

2 结果讨论

2.1 计算参数

在本节中，展示了250Cf 诱发裂变的研究结果，其中裂变碎片的电荷分布显示对称裂变峰占据主导地位。在第一步中，执行大规模形变约束自洽计算以生成(β2，β3)平面中的势能面和单核子波函数。集体变量β2的范围为(-1,6.0)，β3的范围为(0.0,3.12)，计算步长Δβ2=Δβ3=0.04。能量密度泛函采用PC-PK1 参数组，δ- 对力的强度参数采用数值：Vn=349.5MeVfm3以及Vp=330MeVfm3。在轴向形变的谐振子基上求解了单粒子波函数的自洽Dirac 方程。采用Felix2.0程序用于模拟裂变核的时间演化，时间步长δt=5*10-4zs。在剪裂线之外的区域，考虑集体波包逃逸的附加虚吸收势参数为：吸收率r=30*1022s-1，吸收带宽度w=1.5。

图1(β2，β3)平面内250Cf 的自洽RMF+BCS 四极和八极约束形变位能曲面(单位：MeV)

2.2 计算结果

图1 显示了自洽RMF+BCS 四极和八极约束能量位能曲面以及部分点的密度分布。在基态谷内，能量最低点位于(β2，β3)～(0.33，0.00)附近。整个位能曲面主要存在一个非对称裂变谷。在(β2，β3)～(0.50，0.00) 以及(β2，β3)～(1.88，0.30)附近存在势垒。发生裂变时，集体波函数受非对称裂变谷的影响将会向大β3的方向演化，从微观角度预测了250Cf 的裂变以非对称裂变为主。位能曲面左上角同时给出了结团发射的反应道。

图2(β2，β3)平面内集体波函数模方|g|2 在不同时刻时的分布图(含剪裂线)。波函数由剪裂线内动态演化流出剪裂线，流出后认为裂变事件发生

在描述(β2，β3)平面的裂变时，判断裂变是否发生的一个关键因素是裂变区域与裂变后区域的不连续性，这个不连续性在能量、脖子数、密度分布等各个物理量上都有体现。定义脖子数算符：

其中zN为脖子位置，aN取为1。遵循文献[20]的结果，我们拟定＜Q^N＞=3.0 为剪裂线。

类薛定谔方程描述了集体波函数在四八极集体空间中的演化规律，它和初始条件(初态)，边界条件一起决定了最终的产物分布。在本文中，对于剪裂线外部分的边界条件处理采用Felix2.0 程序自带的方式[14]进行外推。而初态则采用boost(冲量近似态)的方式构建。

当初态与边界确定后，则可以根据类薛定谔方程追踪波函数及其模方在集体空间中的演化。

在图2 中，我们展示了不同时刻集体波函数模方|g|2 的分布图(初态能量高于裂变势垒1MeV)。T=0zs 对应波函数演化的初态，其主要分布都存在于基态谷内。随着时间的流逝，可以看见，波函数开始朝着大形变方向演化。由于(β2，β3)～(0.33，0.00)附近的裂变势垒的影响，可以看见图2(b)中波函数分为了二个分支，分别沿着正负β3。这种演化规律一直持续至裂变完成。从图2(c)和图2(d)中可以看见，非对称裂变谷处的波函数流出量依旧占据主导地位。值得注意的是，虽然对称裂变谷能量较高，但由于在集体空间中距离基态的距离较近，依旧存在不可忽视的集体流通量。

另一个值得注意的量是剪裂线上波函数通量的累计值。如图3 所示，随着时间流逝，总通量先维持0 值一小段时间，后快速增加，然后趋于平缓。裂变的前2.5zs，对应着基态谷中的集体波函数受激后越垒的过程，此时并没有集体波函数从剪裂线上流出，但实际上原子核正发生着剧烈的形状演化。之后，累计通量快速上升，代表着集体波函数大量流出，这意味着原子核裂变事件发生的概率快速上升。最后在40zs以后开始趋于平缓。在这之后虽然依旧有波函数流出，但最终的裂变产物分布已经趋于稳定。

图3 剪裂线上波函数通量的累计值，波函数总量归一。随着时间流逝，总通量快速增加后缓慢收敛

最后，我们给出的裂变产物分布(图4)也佐证了此点。需要指出的是，虽然不同时刻从剪裂线上流出的波函数总量的绝对数值存在差异，但经过高斯平滑和归一化后，不同时刻的产物分布在图4 上的绝对面积相同，而只是不同产物的分布占比不同。正如我们预料的一样，由于对称裂变谷的距离过近，裂变的前5 zs 为对称裂变占据主导。之后随着时间演化，对称峰下降，在10 zs 时整个产额分布便已达到平衡，之后虽然依旧有大量波函数流出，但对产额的分布已无影响。

图4 250Cf 的中子发射前裂变碎片分布。分布数据都归一化到200%

3 结论

利用协变密度泛函理论计算了250Cf 的静态位能曲面，通过BCS 近似考虑了核子的对关联，通过微扰推转近似得到了原子核的集体质量张量，通过原子核总能量扣除零点能的方式得到了集体势场，并使用含时生成坐标方法和高斯重叠近似分析了集体空间中的波函数演化。

TDGCM+GOA 计算再现了裂变质量分布的主要特征。我们发现，由于位能曲面的独特拓朴结构，250Cf 的低能诱发裂变呈现四峰结构，且产额分布在10zs 内便已收敛，后续演化只会带来剪裂线上流出的总通量数值增加，对产额分布并无影响。