234U裂变的五维势能曲面计算

朱 鑫，王智明，罗昶恺，倪 磊，钟春来，樊铁栓,*

(1.集美大学 理学院，福建 厦门 361021；2.北京大学 物理学院 重离子研究所 核物理与核技术国家重点实验室，北京 100871)

在处于激发态的原子核通过大形变的集体运动分裂成两个碎片并释放中子的过程中，核势能随形变参量的变化可由一多维空间上的宏观-微观势能函数(通常称为势能曲面)表征。势能曲面反映了核体系能量在裂变过程中不断变化的状态，对于理解裂变机制以及分析裂变产物的性质具有决定性意义，因此是裂变的宏观-微观理论和主要的微观理论模型的重要物理基础[1-3]。微观理论的代表性理论[4-7]包括核密度泛函理论(density functional theory, DFT)和自洽平均场(self-consistent mean field, SCMF)理论，如能量密度泛函(EDF)方法[8]。在宏观-微观模型建立和发展过程中，Brosa等[9-10]建立了随机颈部断裂模型，并首次获得了五维变形空间中的裂变势能面，首次从三维势能曲面上发现存在多裂变通道，主要包括标准道、超长道和超短通道，并在后来出版的著作中得到证实[11-13]。另一方面，Möller等[14-17]利用有限程液滴模型(FRLDM)计算裂变宏观能量，单粒子势能则采用有限程折叠汤川势模型计算，给出质子数从78到125的共1 585个核的裂变势垒计算结果[17]，也证实了在大多数重核中存在对称和不对称裂变的两种裂变模式[18]。近20年来，包括双中心壳模型在内的多种裂变势能曲面计算方法与裂变动力学结合在裂变后现象的定量预测上连续取得了重要突破[19-23]。

宏观-微观模型的一个重要问题是进行合理的形状参数化，所选择的集体坐标系统应能够准确描述从基态到断裂构型的核形状演变过程。由于过多的空间自由度所增加的计算成本将不可避免地抑制该参数化的可访问性，所以现实的形变参数选择应既准确又简单，以便将计算量减少到合理。本文将采用5个形变参量，以计算得到五维势能曲面。

本文利用宏观-微观模型计算234U复合核裂变的五维势能曲面，采用搜索程序在五维势能曲面上获得裂变势垒以及鞍点形状。

1 五维势能曲面的计算

在计算宏观能量时使用考虑了核表面曲率效应的液滴模型的LSD(Lublin-Strasbourg drop)公式[24]，微观能的计算中选择折叠汤川势[25-28]作为单粒子势。微观修正使用Strutinsky壳修正方法[29-30]和SBCS(Strutinsky smoothed BCS)对修正方法[31]。选取Nix连接二次曲面核形状[32]进行从椭球形的基态到哑铃状的断前形状描述，针对每个形状计算宏观能与微观能来得到对应的形变核势能。

1.1 宏观-微观模型

在宏观-微观模型中，作为形变参数函数的总势能E由宏观能Emac和微观修正能Emic相加组成：

E=Emac+Emic

(1)

在略去与形变无关的部分后，采用的LSD模型公式[24]可写为：

Emac=bsurf(1-κsurfI2)A2/3Bsurf+

bcur(1-κcurI2)A1/3Bcur+

(2)

其中：I=|N-Z|/A,为裂变核的同位旋，N、Z和A分别为裂变核的中子数、质子数和质量数；Bsurf、Bcur和Bcoul分别为变形核的表面积、曲率积分和库仑能与球形核对应量的比值；拟合系数分别选为bsurf=16.970 7 MeV,κsurf=2.293 8，bcur=3.860 2 MeV,κcur=-2.476 4,r0=1.217 25 fm。

式(2)中的系数是通过对2 766个原子核结合能的实验数据拟合得到的。虽然这些系数并未对裂变位垒的实验数据进行拟合，但经过验证，LSD公式可很好应用于裂变势垒的计算。

以宏观能计算为基础，考虑核壳层结构和成对效应，本工作采用Strutinsky方法进行微观能量修正。在微观修正能计算中，采用折叠汤川势作为单粒子势[25-28]：

(3)

式中：a=0.80 fm，为折叠程；V0为势阱深度，对中子和质子取值不同。

根据一种近似的变分方法(BCS)可求得因为对效应而附加的能量，从而进行对修正。本工作的计算采用Olofsson等[31]的SBCS公式。

1.2 核形状描述与形变参量

在宏观-微观模型计算中，常用的形状描述方法包括广义椭球[32-34]、卡西尼卵形线[35]和连接二次曲面。连接多个曲面以描述原子核形状的想法最初是由Grammaticos等[36]提出的。后来，Swiatecki意识到伸长变形的原子核可能有一颈部，他在连接曲面中间增加了一圆柱形的颈部[37]，得到了一类似哑铃形的连接曲面，Nix进而将哑铃形颈部用二次曲面代替[38-39]。本文用3个二次曲面光滑连接的方式描述形状演化(图1),当复合核拉长至越过外位垒后,中间的二次曲面由椭球形变成了双曲形,进而出现哑铃状形变：

图1 连接二次曲面示意图Fig.1 Visualization of three quadratic surfaces

(4)

其中：a、c和l分别为长半轴长、短半轴长和球体中心的位置；下标1和2分别表示左侧曲面和右侧曲面，这两个曲面都是椭球体；中间的椭球体由下标3标记，c3为实数时代表球体，c3为虚数时代表双曲面；z1和z2分别为中间球体与两个端部球体的左、右连接点。在实际运算中对此形状进行平移操作，将坐标原点从任意位置移动到了l3=0处。然而，这些参数并不是独立的，通过适当的平滑连接约束只剩下5个独立的参数。

Nix在论述这些形状参数的文章中提出了1组用几何参量表示的对称参量(σ1、σ2、σ3)与非对称参量(α1、α2、α3)，根据约束条件可将几何参量用5个独立参量表示[38]。

为更真实反映物理，进而使用电四极矩Q2、颈部半径d、不对称参数αg、左侧碎片变形参数εf1和右侧碎片变形参数εf2作为1组变形参数描述核形状的演化[15]。Nix的两组形变参量可作为连接这组形变参量与几何参量的纽带。

其中，电四极矩的表达式为：

(5)

将非长度量纲的电荷及其他系数忽略后，本工作采用的具体公式为：

(6)

而颈部半径d有：

d=a3

(7)

另3个变量的定义取自Moller等[15]的工作，表达式分别为：

(8)

(9)

(10)

其中，M1、M2分别为左、右碎片质量。在连接二次曲面参数化中，五维势能面上网格点的密度选择应保持适中，这是由于壳修正方法计算的原子核结合能存在一定波动，当网格足够细时，进一步减小步长不会使结果更准确，反而会浪费大量计算时间。若网格间距过小，关键点(基态、鞍点和断点)和从中提取的路径信息反而不够准确。在这项工作中，本工作经过大量测试与调试后选择的网格步长如下：

(11)

αg=0.02(i-2)i=1,2,…,35

(12)

d=0.08Rcn(16-i)i=1,2,…,15

(13)

εf1=(-0.2,-0.15,-0.1,0.00,0.1,

0.15,0.175,0.2,0.225,0.25,0.275,

0.3,0.35,0.4,0.5)

(14)

εf2=(-0.2,-0.15,-0.1,0.00,

0.1,0.15,0.175,0.2,0.225,0.25,

0.275,0.3,0.35,0.4,0.5)

(15)

其中，Rcn为球形核半径。所以，本工作最后所选择的格点总数为35×50×15×15×15=5 906 250。

2 多模式裂变路径

对于多模式的裂变路径，通常是投影到其中最重要的两个维度上寻找近似的裂变路径。对于复合核234U，将计算得到的五维势能曲面对Q2和αg两个方向进行投影，得到二维立体投影图(图2)与二维地形图(图3)。由Edef(Q2,d,αg,εf1,εf2)对其他3个坐标(d,εf1,εf2)取最小值得到：

图2 复合核234U的势能曲面计算结果在Q2和αg方向的二维投影图Fig.2 Two-dimensional projection of five-dimensional potential energy surface in Q2 and αg directions for 234U

图3 复合核234U的势能曲面计算结果在Q2和αg方向的二维投影地形图Fig.3 Two-dimensional projection contour map of five-dimensional potential energy surface in Q2 and αg directions for 234U

(16)

如图3所示，复合核234U存在两个分离良好的裂变路径，共享相同的内势垒EA，并在第2个势阱处向下偏离。不对称路径延伸至A点结束，而对称路径在S点结束。可发现，不对称参数αg约11.3，根据式(10)，这意味着重碎片的质量数接近140。而对于对称路径，αg约为0。图中，gs表示基态的近似位置。基态重核的形状应是椭球体，所以αg应约为0。第2个势垒即裂变的外势垒在图中标记为EB，其能量值为外势垒的高度。两条裂变路径在大形变时明显分开，快到达断裂点时，非对称裂变路径将穿过一比外势垒EB高度更低的新势垒。但对称裂变路径需穿过的势垒高度高于EB，所以，对称裂变路径的外势垒高度应是需再次穿过的势垒高度。

3 势垒高度

作为决定裂变反应与其他衰变道竞争中胜出概率的一关键物理量，从五维势能曲面计算中得到裂变势垒的性质具有重要的应用意义。在势能曲面计算中，一旦确定了内鞍点和外鞍点的能量，减去基态能量即可得到势垒高度，同时可得到势能曲面上这些鞍点的形变参数，再通过形变参量得到几何参数，从而画出这些位置对应的核形状。然而，在五维势能面上搜索数百万个格点是一项艰巨的任务，因此，使用投影后的二维势能曲面上的鞍点位置作为内、外鞍点的位置。

本工作采用模拟泛洪算法[40]对复合核234U的势能曲面进行搜索，可得到4 134个极小值点,将投影的每个Edef(Q2,αg)的其余3个坐标输出，就得到了Edef(Q2,d,αg,εf1,εf2)。在复合核234U势能曲面的二维投影上寻找附近的点，得到一对极小值点(10,1,6,5,7)与(10,3,6,7,5)，认为这对互为镜像的极小值点对应的核形状是234U裂变核内鞍点的形状。外鞍点的近似形状选择了(21,20,8,4,11)点的形状，另外第2势阱近似形状即为极小值点(17,9,6,4,13)的形状。以上这些形状与基态近似形状均在图4中给出。

图4 复合核234U的近似基态(a)、近似内鞍点(b)、近似第2势阱(c)、近似外鞍点(d)的核形状Fig.4 Approximate nuclear shape in ground states (a), inner saddle point (b), second potential well (c) and outer saddle point (d) of 234U

从这些特殊位置的近似形状图上可看出，连接二次曲面对于大形变的描述良好，也从一个方面显示了本工作势能曲面的计算结果的可靠性。然而，对于基态的位置，连接二次曲面形状描述并不是很好，这是由该形状描述本身的特点决定的。只有当a1=a2=a3以及c1=c2=c3时，3段椭球才能完美地连成1个椭球，给出基态的形状。因为体积固定，所以对应的五维形变参量也是唯一确定的。但这个形状很可能并不在本工作给出的格点中，若想尝试给出每个裂变核的基态形状，就必须对每个核修改五维形变参量的格点步长，这样就难以给出统一的格点步长。

本工作采用二维投影寻找裂变路径上的特征极小值点作为模拟泛洪搜索程序[40]的入水点与出水点，搜索了内、外两个鞍点的能量。考虑到连接二次曲面形状对基态描述的缺陷，并且对于不同形变区域采用不同核形状描述被证明是可行的[41]。因此，将计算得到的内、外两个鞍点的能量减去由广义Lawrence形状描述计算得到的基态能量[41-44]，就得到了内、外势垒的高度。复合核234U基态能量Egs为-3.643 MeV，内鞍点能量Edef,in为0.717 MeV，外鞍点能量Edef,out为0.677 MeV，内势垒高度EA为4.360 MeV，外势垒高度EB为4.320 MeV。

本工作计算的234U裂变势垒高度与实验数据和其他理论结果的比较列于表1。其中，Moller等[15]使用连接二次曲面进行形状描述，而Zhong等[42]采用的是广义Lawrence形状描述。Wang等[44]在广义Lawrence形状描述下采用双中心谐振子基矢完成了哈密顿矩阵计算。本工作的内垒高度与实验数据最接近，但外垒高度最低，未来也可通过将单中心谐振子基矢换成双中心谐振子基矢改进本工作结果。

表1 复合核234U势垒高度的比较Table 1 Comparison of barrier heights for compound nucleus 234U

4 结论

应用宏观-微观方法计算了复合核234U的包含5 906 250个核形状的五维势能曲面。宏观-微观模型由连接二次曲面进行核形状描述，采用LSD公式加微观壳修正的方法计算形变核的势能，单粒子势选用折叠汤川势。在二维势能曲面投影上找出对称裂变和非对称裂变两个裂变通道，并得到各特殊位置对应五维势能曲面上格点的5个形变参量，得到了5 906 250个不同的形变核形状。在五维势能曲面上采用裂变势垒搜索程序，获得了复合核234U的裂变势垒，所得到的势垒高度与其他理论结果和RIPL库数据符合较好。基于该裂变位能曲面的五维朗之万动力学研究裂变后现象的工作正在进行中。