234U裂變的五維勢能曲面計算

朱 鑫，王智明，羅昶愷，倪 磊，鐘春來，樊鐵栓,*

(1.集美大學(xué) 理學(xué)院，福建 廈門 361021；2.北京大學(xué) 物理學(xué)院 重離子研究所 核物理與核技術(shù)國家重點實驗室，北京 100871)

在處于激發(fā)態(tài)的原子核通過大形變的集體運動分裂成兩個碎片并釋放中子的過程中，核勢能隨形變參量的變化可由一多維空間上的宏觀-微觀勢能函數(shù)(通常稱為勢能曲面)表征。勢能曲面反映了核體系能量在裂變過程中不斷變化的狀態(tài)，對于理解裂變機(jī)制以及分析裂變產(chǎn)物的性質(zhì)具有決定性意義，因此是裂變的宏觀-微觀理論和主要的微觀理論模型的重要物理基礎(chǔ)[1-3]。微觀理論的代表性理論[4-7]包括核密度泛函理論(density functional theory, DFT)和自洽平均場(self-consistent mean field, SCMF)理論，如能量密度泛函(EDF)方法[8]。在宏觀-微觀模型建立和發(fā)展過程中，Brosa等[9-10]建立了隨機(jī)頸部斷裂模型，并首次獲得了五維變形空間中的裂變勢能面，首次從三維勢能曲面上發(fā)現(xiàn)存在多裂變通道，主要包括標(biāo)準(zhǔn)道、超長道和超短通道，并在后來出版的著作中得到證實[11-13]。另一方面，M?ller等[14-17]利用有限程液滴模型(FRLDM)計算裂變宏觀能量，單粒子勢能則采用有限程折疊湯川勢模型計算，給出質(zhì)子數(shù)從78到125的共1 585個核的裂變勢壘計算結(jié)果[17]，也證實了在大多數(shù)重核中存在對稱和不對稱裂變的兩種裂變模式[18]。近20年來，包括雙中心殼模型在內(nèi)的多種裂變勢能曲面計算方法與裂變動力學(xué)結(jié)合在裂變后現(xiàn)象的定量預(yù)測上連續(xù)取得了重要突破[19-23]。

宏觀-微觀模型的一個重要問題是進(jìn)行合理的形狀參數(shù)化，所選擇的集體坐標(biāo)系統(tǒng)應(yīng)能夠準(zhǔn)確描述從基態(tài)到斷裂構(gòu)型的核形狀演變過程。由于過多的空間自由度所增加的計算成本將不可避免地抑制該參數(shù)化的可訪問性，所以現(xiàn)實的形變參數(shù)選擇應(yīng)既準(zhǔn)確又簡單，以便將計算量減少到合理。本文將采用5個形變參量，以計算得到五維勢能曲面。

本文利用宏觀-微觀模型計算234U復(fù)合核裂變的五維勢能曲面，采用搜索程序在五維勢能曲面上獲得裂變勢壘以及鞍點形狀。

1 五維勢能曲面的計算

在計算宏觀能量時使用考慮了核表面曲率效應(yīng)的液滴模型的LSD(Lublin-Strasbourg drop)公式[24]，微觀能的計算中選擇折疊湯川勢[25-28]作為單粒子勢。微觀修正使用Strutinsky殼修正方法[29-30]和SBCS(Strutinsky smoothed BCS)對修正方法[31]。選取Nix連接二次曲面核形狀[32]進(jìn)行從橢球形的基態(tài)到啞鈴狀的斷前形狀描述，針對每個形狀計算宏觀能與微觀能來得到對應(yīng)的形變核勢能。

1.1 宏觀-微觀模型

在宏觀-微觀模型中，作為形變參數(shù)函數(shù)的總勢能E由宏觀能Emac和微觀修正能Emic相加組成：

E=Emac+Emic

(1)

在略去與形變無關(guān)的部分后，采用的LSD模型公式[24]可寫為：

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

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

(2)

其中：I=|N-Z|/A,為裂變核的同位旋，N、Z和A分別為裂變核的中子數(shù)、質(zhì)子數(shù)和質(zhì)量數(shù)；Bsurf、Bcur和Bcoul分別為變形核的表面積、曲率積分和庫侖能與球形核對應(yīng)量的比值；擬合系數(shù)分別選為bsurf=16.970 7 MeV,κsurf=2.293 8，bcur=3.860 2 MeV,κcur=-2.476 4,r0=1.217 25 fm。

式(2)中的系數(shù)是通過對2 766個原子核結(jié)合能的實驗數(shù)據(jù)擬合得到的。雖然這些系數(shù)并未對裂變位壘的實驗數(shù)據(jù)進(jìn)行擬合，但經(jīng)過驗證，LSD公式可很好應(yīng)用于裂變勢壘的計算。

以宏觀能計算為基礎(chǔ)，考慮核殼層結(jié)構(gòu)和成對效應(yīng)，本工作采用Strutinsky方法進(jìn)行微觀能量修正。在微觀修正能計算中，采用折疊湯川勢作為單粒子勢[25-28]：

(3)

式中：a=0.80 fm，為折疊程；V0為勢阱深度，對中子和質(zhì)子取值不同。

根據(jù)一種近似的變分方法(BCS)可求得因為對效應(yīng)而附加的能量，從而進(jìn)行對修正。本工作的計算采用Olofsson等[31]的SBCS公式。

1.2 核形狀描述與形變參量

在宏觀-微觀模型計算中，常用的形狀描述方法包括廣義橢球[32-34]、卡西尼卵形線[35]和連接二次曲面。連接多個曲面以描述原子核形狀的想法最初是由Grammaticos等[36]提出的。后來，Swiatecki意識到伸長變形的原子核可能有一頸部，他在連接曲面中間增加了一圓柱形的頸部[37]，得到了一類似啞鈴形的連接曲面，Nix進(jìn)而將啞鈴形頸部用二次曲面代替[38-39]。本文用3個二次曲面光滑連接的方式描述形狀演化(圖1),當(dāng)復(fù)合核拉長至越過外位壘后,中間的二次曲面由橢球形變成了雙曲形,進(jìn)而出現(xiàn)啞鈴狀形變：

圖1 連接二次曲面示意圖Fig.1 Visualization of three quadratic surfaces

(4)

其中：a、c和l分別為長半軸長、短半軸長和球體中心的位置；下標(biāo)1和2分別表示左側(cè)曲面和右側(cè)曲面，這兩個曲面都是橢球體；中間的橢球體由下標(biāo)3標(biāo)記，c3為實數(shù)時代表球體，c3為虛數(shù)時代表雙曲面；z1和z2分別為中間球體與兩個端部球體的左、右連接點。在實際運算中對此形狀進(jìn)行平移操作，將坐標(biāo)原點從任意位置移動到了l3=0處。然而，這些參數(shù)并不是獨立的，通過適當(dāng)?shù)钠交B接約束只剩下5個獨立的參數(shù)。

Nix在論述這些形狀參數(shù)的文章中提出了1組用幾何參量表示的對稱參量(σ1、σ2、σ3)與非對稱參量(α1、α2、α3)，根據(jù)約束條件可將幾何參量用5個獨立參量表示[38]。

為更真實反映物理，進(jìn)而使用電四極矩Q2、頸部半徑d、不對稱參數(shù)αg、左側(cè)碎片變形參數(shù)εf1和右側(cè)碎片變形參數(shù)εf2作為1組變形參數(shù)描述核形狀的演化[15]。Nix的兩組形變參量可作為連接這組形變參量與幾何參量的紐帶。

其中，電四極矩的表達(dá)式為：

(5)

將非長度量綱的電荷及其他系數(shù)忽略后，本工作采用的具體公式為：

(6)

而頸部半徑d有：

d=a3

(7)

另3個變量的定義取自Moller等[15]的工作，表達(dá)式分別為：

(8)

(9)

(10)

其中，M1、M2分別為左、右碎片質(zhì)量。在連接二次曲面參數(shù)化中，五維勢能面上網(wǎng)格點的密度選擇應(yīng)保持適中，這是由于殼修正方法計算的原子核結(jié)合能存在一定波動，當(dāng)網(wǎng)格足夠細(xì)時，進(jìn)一步減小步長不會使結(jié)果更準(zhǔn)確，反而會浪費大量計算時間。若網(wǎng)格間距過小，關(guān)鍵點(基態(tài)、鞍點和斷點)和從中提取的路徑信息反而不夠準(zhǔn)確。在這項工作中，本工作經(jīng)過大量測試與調(diào)試后選擇的網(wǎng)格步長如下：

(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為球形核半徑。所以，本工作最后所選擇的格點總數(shù)為35×50×15×15×15=5 906 250。

2 多模式裂變路徑

對于多模式的裂變路徑，通常是投影到其中最重要的兩個維度上尋找近似的裂變路徑。對于復(fù)合核234U，將計算得到的五維勢能曲面對Q2和αg兩個方向進(jìn)行投影，得到二維立體投影圖(圖2)與二維地形圖(圖3)。由Edef(Q2,d,αg,εf1,εf2)對其他3個坐標(biāo)(d,εf1,εf2)取最小值得到：

圖2 復(fù)合核234U的勢能曲面計算結(jié)果在Q2和αg方向的二維投影圖Fig.2 Two-dimensional projection of five-dimensional potential energy surface in Q2 and αg directions for 234U

圖3 復(fù)合核234U的勢能曲面計算結(jié)果在Q2和αg方向的二維投影地形圖Fig.3 Two-dimensional projection contour map of five-dimensional potential energy surface in Q2 and αg directions for 234U

(16)

如圖3所示，復(fù)合核234U存在兩個分離良好的裂變路徑，共享相同的內(nèi)勢壘EA，并在第2個勢阱處向下偏離。不對稱路徑延伸至A點結(jié)束，而對稱路徑在S點結(jié)束?？砂l(fā)現(xiàn)，不對稱參數(shù)αg約11.3，根據(jù)式(10)，這意味著重碎片的質(zhì)量數(shù)接近140。而對于對稱路徑，αg約為0。圖中，gs表示基態(tài)的近似位置。基態(tài)重核的形狀應(yīng)是橢球體，所以αg應(yīng)約為0。第2個勢壘即裂變的外勢壘在圖中標(biāo)記為EB，其能量值為外勢壘的高度。兩條裂變路徑在大形變時明顯分開，快到達(dá)斷裂點時，非對稱裂變路徑將穿過一比外勢壘EB高度更低的新勢壘。但對稱裂變路徑需穿過的勢壘高度高于EB，所以，對稱裂變路徑的外勢壘高度應(yīng)是需再次穿過的勢壘高度。

3 勢壘高度

作為決定裂變反應(yīng)與其他衰變道競爭中勝出概率的一關(guān)鍵物理量，從五維勢能曲面計算中得到裂變勢壘的性質(zhì)具有重要的應(yīng)用意義。在勢能曲面計算中，一旦確定了內(nèi)鞍點和外鞍點的能量，減去基態(tài)能量即可得到勢壘高度，同時可得到勢能曲面上這些鞍點的形變參數(shù)，再通過形變參量得到幾何參數(shù)，從而畫出這些位置對應(yīng)的核形狀。然而，在五維勢能面上搜索數(shù)百萬個格點是一項艱巨的任務(wù)，因此，使用投影后的二維勢能曲面上的鞍點位置作為內(nèi)、外鞍點的位置。

本工作采用模擬泛洪算法[40]對復(fù)合核234U的勢能曲面進(jìn)行搜索，可得到4 134個極小值點,將投影的每個Edef(Q2,αg)的其余3個坐標(biāo)輸出，就得到了Edef(Q2,d,αg,εf1,εf2)。在復(fù)合核234U勢能曲面的二維投影上尋找附近的點，得到一對極小值點(10,1,6,5,7)與(10,3,6,7,5)，認(rèn)為這對互為鏡像的極小值點對應(yīng)的核形狀是234U裂變核內(nèi)鞍點的形狀。外鞍點的近似形狀選擇了(21,20,8,4,11)點的形狀，另外第2勢阱近似形狀即為極小值點(17,9,6,4,13)的形狀。以上這些形狀與基態(tài)近似形狀均在圖4中給出。

圖4 復(fù)合核234U的近似基態(tài)(a)、近似內(nèi)鞍點(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

從這些特殊位置的近似形狀圖上可看出，連接二次曲面對于大形變的描述良好，也從一個方面顯示了本工作勢能曲面的計算結(jié)果的可靠性。然而，對于基態(tài)的位置，連接二次曲面形狀描述并不是很好，這是由該形狀描述本身的特點決定的。只有當(dāng)a1=a2=a3以及c1=c2=c3時，3段橢球才能完美地連成1個橢球，給出基態(tài)的形狀。因為體積固定，所以對應(yīng)的五維形變參量也是唯一確定的。但這個形狀很可能并不在本工作給出的格點中，若想嘗試給出每個裂變核的基態(tài)形狀，就必須對每個核修改五維形變參量的格點步長，這樣就難以給出統(tǒng)一的格點步長。

本工作采用二維投影尋找裂變路徑上的特征極小值點作為模擬泛洪搜索程序[40]的入水點與出水點，搜索了內(nèi)、外兩個鞍點的能量?？紤]到連接二次曲面形狀對基態(tài)描述的缺陷，并且對于不同形變區(qū)域采用不同核形狀描述被證明是可行的[41]。因此，將計算得到的內(nèi)、外兩個鞍點的能量減去由廣義Lawrence形狀描述計算得到的基態(tài)能量[41-44]，就得到了內(nèi)、外勢壘的高度。復(fù)合核234U基態(tài)能量Egs為-3.643 MeV，內(nèi)鞍點能量Edef,in為0.717 MeV，外鞍點能量Edef,out為0.677 MeV，內(nèi)勢壘高度EA為4.360 MeV，外勢壘高度EB為4.320 MeV。

本工作計算的234U裂變勢壘高度與實驗數(shù)據(jù)和其他理論結(jié)果的比較列于表1。其中，Moller等[15]使用連接二次曲面進(jìn)行形狀描述，而Zhong等[42]采用的是廣義Lawrence形狀描述。Wang等[44]在廣義Lawrence形狀描述下采用雙中心諧振子基矢完成了哈密頓矩陣計算。本工作的內(nèi)壘高度與實驗數(shù)據(jù)最接近，但外壘高度最低，未來也可通過將單中心諧振子基矢換成雙中心諧振子基矢改進(jìn)本工作結(jié)果。

表1 復(fù)合核234U勢壘高度的比較Table 1 Comparison of barrier heights for compound nucleus 234U

4 結(jié)論

應(yīng)用宏觀-微觀方法計算了復(fù)合核234U的包含5 906 250個核形狀的五維勢能曲面。宏觀-微觀模型由連接二次曲面進(jìn)行核形狀描述，采用LSD公式加微觀殼修正的方法計算形變核的勢能，單粒子勢選用折疊湯川勢。在二維勢能曲面投影上找出對稱裂變和非對稱裂變兩個裂變通道，并得到各特殊位置對應(yīng)五維勢能曲面上格點的5個形變參量，得到了5 906 250個不同的形變核形狀。在五維勢能曲面上采用裂變勢壘搜索程序，獲得了復(fù)合核234U的裂變勢壘，所得到的勢壘高度與其他理論結(jié)果和RIPL庫數(shù)據(jù)符合較好?；谠摿炎兾荒芮娴奈寰S朗之萬動力學(xué)研究裂變后現(xiàn)象的工作正在進(jìn)行中。