基于協(xié)變密度泛函理論的258Fm誘發(fā)裂變動力學研究

陳盛遠，李澤宇，陳永靜，李志攀,*

(1.西南大學 物理科學與技術學院，重慶 400715；2.中國原子能科學研究院 核數(shù)據重點實驗室，中國核數(shù)據中心，北京 102413)

原子核裂變在能源、國防、工農醫(yī)等領域中扮演著至關重要的角色，一直是核物理以及工程物理中的重要研究對象。在基礎研究方面，裂變是決定超重元素穩(wěn)定性的主要因素之一[1]；決定著星系演化中重元素核合成r-過程的終點以及所合成核素的豐度[2]；也是產生短壽命奇特原子核的重要途徑。因此，原子核裂變研究是當今核物理的重要前沿課題[3-4]。

重核裂變是低能核物理中最復雜的問題之一[3-5]。由于重核中通常包含著數(shù)以百計的核子通過多體相互作用耦合在一起，使用第一性原理計算(abinitiocalculation)或相互作用殼模型(interacting shell-model)研究如此復雜的多體問題尚不太現(xiàn)實。因此，當前對于重核裂變的微觀研究通?；谠雍四芰棵芏确汉碚?density functional theory, DFT)。原子核DFT同時考慮了Pauli不相容原理以及核子-核子間相互作用，且自洽地包含了殼效應和量子多體效應，能給出結構更復雜的裂變位壘和斷點性質。近年來，基于DFT對重核裂變位能曲面及裂變位壘開展了一系列的研究工作[6-19]，深化了人們對裂變機制的微觀認知。

然而，為研究裂變動力學性質并定量描述各類觀測量，必須在DFT中引入含時動力學演化，即發(fā)展時間依賴的密度泛函理論(TDDFT)[20-21]。一般情況下，TDDFT計算只能給出原子核裂變的最可幾路徑，并不能描述裂變全過程和產物分布。最近，人們基于TDDFT開展了一些嘗試以給出裂變產物分布，并取得了卓有成效的進展[22]。

含時生成坐標方法(time-dependent generator coordinate method, TDGCM)將原子核波函數(shù)描述為集體形變空間中內稟態(tài)的疊加，能自洽包含裂變過程中的量子漲落效應，是一種完全量子化的方法?；贒FT，通過引入絕熱近似和高斯重疊近似(Gaussian overlap approximation, GOA)[4-23]，構建了可用于實際計算的TDGCM+GOA框架并編寫了相應的計算程序[24-28]。在該框架中，含時Hill-Wheeler方程退化為集體空間中局域的含時類薛定諤方程，其動力學性質完全由DFT計算得到的多維位能曲面和質量參數(shù)等決定，求解該方程即可得到裂變碎片分布等觀測量。該方法最早可追溯到1991年，Berger等提出了時間依賴的集體模型處理原子核裂變，用TDGCM+GOA計算了238U核的裂變碎片動能以及質量分布[24]。隨后，基于非相對論Skyrme或Gogny泛函，采用TDGCM+GOA研究了一些錒系核的非對稱裂變[25-33]。如Regnier等同時采用Skyrme和Gogny泛函很好地再現(xiàn)了中子誘發(fā)239Pu裂變產物分布，并探討了不同的裂變初態(tài)、質量參數(shù)等對計算結果的影響[30]。最近，粒子數(shù)投影、角動量投影、熱漲落等被相繼引入TDGCM+GOA框架[33-35]，進一步改善了模型對實驗結果的定量描述能力并拓展了使用范圍。

協(xié)變密度泛函理論(covariant DFT, CDFT)在研究原子核結構性質上取得了很大成功[36]，其包含諸多優(yōu)點，如標量、矢量密度相互競爭給出新的飽和機制以及自然給出自旋-軌道耦合勢等。近年來，CDFT被廣泛應用于原子核裂變的研究中，很好地再現(xiàn)了重核的裂變位壘結構及自發(fā)裂變壽命等[10,11,17,37-39]。2017年，基于CDFT構建了TDGCM+GOA并詳細研究了226Th光致誘發(fā)裂變動力學性質，定性再現(xiàn)了裂變產物分布三峰結構，并探討了對關聯(lián)強弱對產物分布的影響[39]。隨后，模型被推廣至有限溫度并系統(tǒng)研究了錒系區(qū)原子核的非對稱裂變[34]。

本文將采用基于CDFT的TDGCM+GOA詳細分析258Fm對稱裂變，包括靜態(tài)位能曲面、剪裂線及其組態(tài)的核子密度分布、總動能分布、集體空間波函數(shù)演化及裂變碎片質量分布等，并探討剪裂線的不同判據對總動能及碎片質量分布的影響。

1 基于協(xié)變密度泛函的含時生成坐標方法

低能核裂變可近似描述為1個由少量集體自由度決定的緩慢絕熱過程，本文將采用TDGCM+GOA模型開展研究。在絕熱近似下，原子核的內稟核子自由度與集體形變自由度退耦合，因此計算可分為相對獨立的兩步。第1步，采用約束的CDFT計算原子核在不同形變下的內稟組態(tài)，以得到原子核的位能曲面、集體質量、剪裂線、碎片核子數(shù)等；第2步，基于CDFT計算得到的集體參量構建四極-八極形變空間的含時類薛定諤方程，求解即可得到集體波函數(shù)的動力學演化及裂變產物分布等。關于模型的詳細介紹參見文獻[40]。

1.1 約束CDFT求解集體參量

為研究裂變，需對CDFT進行如下約束計算：

(1)

(2)

(3)

式中，R0=r0A1/3，為原子核半徑，取r0=1.2。

由約束CDFT出發(fā)，自洽求解單核子Dirac方程并結合對關聯(lián)BCS方程便可得到不同形變(β2,β3)下原子核總能量、單核子波函數(shù)及占據幾率等，進而可計算原子核微觀集體質量Bkl(k,l=2,3表示原子核在四極形變和八極形變方向上的集體慣量)和集體勢場Vcoll。通過微擾推轉近似[41-42]，集體質量可表示為：

(4)

M(n),kl(β2,β3)=

(5)

Vcoll(β2,β3)=Etot(β2,β3)-ΔVvib(β2,β3)-

ΔVrot(β2,β3)=Etot(β2,β3)-

(6)

式中，I為轉動慣量，由Inglis-Belyaev公式[45-46]計算得到。

1.2 含時類薛定諤方程

以原子核形變參數(shù)β2和β3為生成坐標，采用高斯重疊近似，從含時Hill-Wheeler方程出發(fā)可推導出集體空間中局域的含時類薛定諤方程：

(7)

式中，g(β2,β3,t)為在(β2,β3)形變空間中的集體波函數(shù)，集體參量Bkl和Vcoll皆由CDFT自洽計算得到。求解該方程即可得集體波函數(shù)隨時間的演化，進而可計算概率流密度：

(8)

為計算碎片分布還需定義剪裂線，即原子核拉伸到剛好要斷裂時所對應的形變組態(tài)的集合。首先定義脖子處核子數(shù)算符：

(9)

(10)

即可得到該線元對應斷裂組態(tài)發(fā)生的概率，考慮整個剪裂線即得到最終產物分布。

2 結果討論

本節(jié)詳細分析258Fm對稱裂變的理論計算結果，數(shù)值細節(jié)如下：CDFT計算采用相對論PC-PK1泛函[47]，對關聯(lián)由BCS方法描述，其中δ-對力的強度參數(shù)通過再現(xiàn)由五點公式[48]計算得到的奇偶質量差給出：Vn=351 MeV·fm3，Vp=366 MeV·fm3。形變約束計算的范圍為β2∈[-0.96,6.0]，β3∈[0.00,3.12]，步長為Δβ2=Δβ3=0.04。采用Felix2.0程序[28]模擬核裂變隨時間的演化過程，時間步長δt=5×10-4zs。在剪裂線外的區(qū)域，考慮集體波包逃逸的附加虛吸收勢參數(shù)為：吸收率r=30×1022s-1，吸收帶寬度w=1.5。

2.1 位能曲面與剪裂線

核密度分布圖由下至上對應的(β2,β3)分別為(2.36, 0.00)、(3.24, 0.76)、(4.28, 1.60)和(4.36, 2.24)圖1 約束CDFT計算得到的258Fm裂變位能曲面在(β2,β3)平面上的分布Fig.1 Potential energy surface of 258Fm in (β2,β3) plane calculated by constrained CDFT

圖1給出了CDFT約束計算得到的258Fm裂變位能曲面在(β2,β3)空間的分布，以及部分形變組態(tài)(β2,β3)=(2.36,0.00),(3.24,0.76),(4.28,1.60),(4.36,2.24)對應的核子密度分布。在第1個勢阱內，能量最低點位于(β2,β3)約(0.28,0.00)附近。在(β2,β3)約(0.68,0.00)處存在裂變位壘，壘高7.3 MeV。在大形變處位能曲面處展現(xiàn)出兩個裂變谷，分別為(β2,β3)約(1.50,0.00)至(3.00,0.00)的對稱裂變谷，以及(β2,β3)約(2.50,0.50)至(4.40,2.24)的非對稱裂變谷。發(fā)生對稱裂變時，由于對稱裂變產物129Sn為幻數(shù)核，子核更傾向于近球形而不會具有較大形變。因此，對稱裂變谷對應的β2較小，僅約1.5，與Gogny D1S有效相互作用的計算結果一致[49]。

由于大形變組態(tài)的復雜性，位能曲面中不同區(qū)域的原子核組態(tài)的脖子處核子數(shù)變化趨勢并不一致，導致不同Qn判據給出的剪裂線存在差異，進而影響裂變總動能及碎片分布。圖2a給出了Qn=4,3,2,1得到的剪裂線在(β2,β3)平面的分布。在非對稱裂變谷附近，Qn=4,3,2的剪裂線明顯分開，而Qn=2,1的剪裂線則近乎重合。這是由于在非對稱裂變谷末端(圖2c)，隨β2增大，脖子處核子數(shù)會平滑減少至2.7左右，再跳變至0.1以下。沿對稱裂變谷(圖2b)，脖子處核子數(shù)由3.2跳變?yōu)?.6后再相對緩慢降至0.1以下，導致Qn=4,3,2時的剪裂線對應的β2取值接近2.38，而Qn=1時的剪裂線對應的β2取值為2.58。雖然剪裂線的不同對集體波函數(shù)的演化影響甚微，但仍可能帶來裂變碎片總動能、碎片產物分布的差異。

圖2 258Fm剪裂線在(β2,β3)平面的分布(a)和β3為0.00(對稱裂變谷)(b)以及2.16(非對稱裂變谷)(c)時脖子處核子數(shù)隨β2的變化關系Fig.2 Scission lines for 258Fm in (β2,β3) plane (a) and evolution of number of nucleons in neck as functions of β2 for β3=0.00 (b) and 2.16 (c)

2.2 裂變碎片總動能

根據裂變剪裂線的定義，可采用庫倫公式估算裂變碎片總動能(TKE)隨核子數(shù)的分布：

(11)

式中：e為電荷單位；dch為兩碎片電荷中心之間的距離；ZH與ZL分別為重、輕碎片的質子數(shù)。

圖3展示了不同Qn判據給出的258Fm初級碎片總動能隨碎片質量數(shù)的分布，作為比較，圖中展示了Faust的計算結果(黑色方塊)[50]，由于目前模型的局限性，暫時無法給出文獻[51]中碎片總動能的概率分布?？倓幽芊植汲拭黠@的單峰結構，這是由于對稱裂變發(fā)生時，兩個碎片距離較近，僅14.5 fm(作為比較，同為錒系核的226Th發(fā)生對稱裂變時，碎片距離為19.86 fm，Qn取3.0)。隨著碎片質量數(shù)偏離平均分裂，碎片之間的相對距離變長(β2增大，見圖2a)，總動能降低。不同Qn對應的總動能分布峰值接近。但Qn較大時的TKE分布更離散，這是由于非對稱裂變時，越大的Qn取值對應著越小的β2形變以及更近的碎片距離。值得注意的是，碎片質量數(shù)122和135附近存在兩組明顯不同的總動能。其原因是(β2,β3)約(3.24,0.76)以及(β2,β3)約(4.28,1.60)附近區(qū)域給出了相似的輕重碎片質量數(shù)(密度分布見圖1)，但斷裂時的碎片距離不同。能量密度泛函在計算裂變位能曲面時，由于剪裂線附近組態(tài)復雜，存在組態(tài)跳變情況，導致沿剪裂線的碎片質量并不連續(xù)，因此，計算得到的總動能存在不連續(xù)的情況，如圖3中給出沿剪裂線重碎片隨八極形變變化的插圖所示。

圖3 不同Qn判據給出的258Fm初級碎片總動能隨碎片質量數(shù)的分布Fig.3 Total kinetic energy distributions of fission fragment mass of 258Fm for different definitions of scission line

2.3 集體空間中波函數(shù)的演化

類薛定諤方程描述了集體波函數(shù)在形變空間中的演化規(guī)律，它和初始條件(初態(tài))、邊界條件一起決定了最終的碎片分布。在本文中，初態(tài)的構建方式與文獻[39]相同，激發(fā)能為8.3 MeV(高于裂變勢壘1 MeV)。剪裂線外部區(qū)域采用能量線性下降的方式進行外推，即以剪裂點為初始點，β2每增大0.04，能量降低1 MeV。圖4給出了不同時刻集體概率密度|g|2在(β2,β3)平面上的分布，紅色線條為Qn=3.0時定義的剪裂線。圖4a對應波函數(shù)演化的初態(tài)，其主要分布于第1個勢阱內。隨著時間推移，波函數(shù)開始朝著大形變方向演化(圖4b)。由圖4c、d可看出，集體波函數(shù)主要于(β2,β3)約(2.36,0.00)附近的對稱裂變道流出，從微觀演化的角度說明了258Fm主要發(fā)生對稱裂變。小部分集體波函數(shù)從(β2,β3)約(4.36,2.24)附近的非對稱裂變道流出。需注意的是，(β2,β3)約(1.40,2.00)附近也有極小部分集體波函數(shù)流出，這是因為位能曲面中同時包含了結團發(fā)射(52Ca)的反應道，但此處流出的波函數(shù)通量相比裂變道可忽略。

2.4 裂變碎片分布

在剪裂線上累計集體流的通量，通過高斯平滑[28]后即可得到裂變碎片分布如圖5所示，藍色實線、紅色短劃線、綠色點劃線和紫色雙點劃線分別對應于剪裂線判據采用Qn=4,3,2,1的產物分布，黑色實心方塊為Flynn于1975年測得的熱中子誘發(fā)裂變的發(fā)射中子后碎片分布數(shù)據[52]，黑色空心圓形和灰色實心方塊分別為Hoffman于1980年和Hulet于1989年測得258Fm自發(fā)裂變的發(fā)射中子后碎片分布數(shù)據[51，53]。整體而言，258Fm以對稱裂變?yōu)橹?，在A約108、150附近理論計算給出矮峰。隨著剪裂線Qn取值由4降至1，對稱裂變峰值由9.88%提升至10.28%，非對稱裂變峰降低并向外偏移。當Qn取4時，非對稱裂變峰在核子數(shù)為110以及148處達到峰值，Qn取1.0時則在108以及150處達到峰值，非對稱裂變峰位移動2個核子。

圖5 TDGCM+GOA計算得到的258Fm誘發(fā)裂變(中子發(fā)射前)碎片質量數(shù)分布并與實驗結果(中子發(fā)射后)[51-53]進行比較Fig.5 Preneutron-emission mass yield for induced fission of 258Fm calculated by TDGCM+GOA in comparison with available data (postneutron-emission)[51-53]

值得注意的是，計算給出的碎片分布峰值介于誘發(fā)裂變實驗數(shù)據與自發(fā)裂變實驗數(shù)據之間，該差異可能與初態(tài)激發(fā)能有關。針對不同激發(fā)能計算得到的碎片分布如圖6所示(Qn=4)。隨著激發(fā)能從8.3 MeV增加至17.3 MeV，集體波函數(shù)在集體空間中的演化變得更加彌散，導致對稱裂變峰值由9.88%降至8.55%。這與實驗數(shù)據中自發(fā)裂變峰值較高，誘發(fā)裂變分布更平緩的表現(xiàn)一致。Regnier等基于Gogny D1S參數(shù)組也研究了Fm同位素鏈裂變動力學[32]，同樣給出了258Fm裂變碎片峰值隨激發(fā)能升高而降低的演化趨勢。與本文相比，文獻[32]中的碎片分布峰值更低，改變激發(fā)能帶來的碎片分布變化更劇烈。

圖6 258Fm誘發(fā)裂變(中子發(fā)射前)碎片質量數(shù)分布隨激發(fā)能的改變并與實驗結果(中子發(fā)射后)[51-53]進行比較Fig.6 Preneutron-emission mass yields for induced fission of 258Fm with different excitation energy in comparison with available data (postneutron-emission)[51-53]

3 總結與展望

基于協(xié)變能量密度泛函PC-PK1計算了258Fm在四極形變與八極形變自由度下的裂變位能曲面，在(β2,β3)約(0.68,0.00)處存在裂變位壘，壘高7.3 MeV，壘外存在路徑較短的對稱裂變谷和路徑較長的非對稱裂變谷。由脖子處粒子數(shù)Qn定義的剪裂線組態(tài)出發(fā)計算得到了裂變碎片總動能分布，其呈現(xiàn)明顯的單峰結構，且隨Qn增大總動能分布變寬。進一步采用TDGCM+GOA方法計算了258Fm低能誘發(fā)裂變時集體波函數(shù)在四極-八極形變空間中的演化，分析了剪裂線和激發(fā)能對裂變碎片質量分布的影響。研究表明，258Fm低能誘發(fā)裂變的碎片質量分布同樣呈現(xiàn)單峰結構，隨著剪裂線Qn取值降低，對稱裂變峰值由9.88%增至10.28%，非對稱裂變峰降低并偏移。同時碎片質量分布受激發(fā)能影響，對稱裂變峰隨激發(fā)能增加而降低，與實驗數(shù)據給出的趨勢一致。

雖然剪裂線帶來的影響并不會導致碎片分布結構的實質性改變，但如何減小甚至消除其帶來的不確定性，特別是對稱峰的高低以及非對稱峰位的改變。一個可能的解決方案是在計算中添加對脖子處核子數(shù)的約束，實現(xiàn)四極、八極以及脖子處核子數(shù)的(β2，β3，Qn)三維約束計算，此外，多維約束的引入有望實現(xiàn)總動能的概率分布。Han的研究表明[54]，在已有四極矩與八極矩約束計算中考慮脖子處核子數(shù)的約束后，得到的多維集體空間中相鄰形變點的密度分布將連續(xù)變化。在新的多維約束空間中，或許可找到更合適的剪裂線定義，從而得到穩(wěn)定且與實驗更符合的結果。