基于協變能量密度泛函的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 內便已收斂，后續演化只會帶來剪裂線上流出的總通量數值增加，對產額分布并無影響。