陳盛遠(yuǎn)，李澤宇，陳永靜，李志攀,*

(1.西南大學(xué) 物理科學(xué)與技術(shù)學(xué)院，重慶 400715；2.中國(guó)原子能科學(xué)研究院 核數(shù)據(jù)重點(diǎn)實(shí)驗(yàn)室，中國(guó)核數(shù)據(jù)中心，北京 102413)

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

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

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

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

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

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

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

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

1.1 約束CDFT求解集體參量

為研究裂變，需對(duì)CDFT進(jìn)行如下約束計(jì)算：

(1)

(2)

(3)

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

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

(4)

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

(5)

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

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

(6)

式中，I為轉(zhuǎn)動(dòng)慣量，由Inglis-Belyaev公式[45-46]計(jì)算得到。

1.2 含時(shí)類薛定諤方程

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

(7)

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

(8)

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

(9)

(10)

即可得到該線元對(duì)應(yīng)斷裂組態(tài)發(fā)生的概率，考慮整個(gè)剪裂線即得到最終產(chǎn)物分布。

2 結(jié)果討論

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

2.1 位能曲面與剪裂線

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

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

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

圖2 258Fm剪裂線在(β2,β3)平面的分布(a)和β3為0.00(對(duì)稱裂變谷)(b)以及2.16(非對(duì)稱裂變谷)(c)時(shí)脖子處核子數(shù)隨β2的變化關(guān)系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 裂變碎片總動(dòng)能

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

(11)

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

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

圖3 不同Qn判據(jù)給出的258Fm初級(jí)碎片總動(dòng)能隨碎片質(zhì)量數(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)的構(gòu)建方式與文獻(xiàn)[39]相同，激發(fā)能為8.3 MeV(高于裂變勢(shì)壘1 MeV)。剪裂線外部區(qū)域采用能量線性下降的方式進(jìn)行外推，即以剪裂點(diǎn)為初始點(diǎn)，β2每增大0.04，能量降低1 MeV。圖4給出了不同時(shí)刻集體概率密度|g|2在(β2,β3)平面上的分布，紅色線條為Qn=3.0時(shí)定義的剪裂線。圖4a對(duì)應(yīng)波函數(shù)演化的初態(tài)，其主要分布于第1個(gè)勢(shì)阱內(nèi)。隨著時(shí)間推移，波函數(shù)開始朝著大形變方向演化(圖4b)。由圖4c、d可看出，集體波函數(shù)主要于(β2,β3)約(2.36,0.00)附近的對(duì)稱裂變道流出，從微觀演化的角度說明了258Fm主要發(fā)生對(duì)稱裂變。小部分集體波函數(shù)從(β2,β3)約(4.36,2.24)附近的非對(duì)稱裂變道流出。需注意的是，(β2,β3)約(1.40,2.00)附近也有極小部分集體波函數(shù)流出，這是因?yàn)槲荒芮嬷型瑫r(shí)包含了結(jié)團(tuán)發(fā)射(52Ca)的反應(yīng)道，但此處流出的波函數(shù)通量相比裂變道可忽略。

2.4 裂變碎片分布

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

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

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

圖6 258Fm誘發(fā)裂變(中子發(fā)射前)碎片質(zhì)量數(shù)分布隨激發(fā)能的改變并與實(shí)驗(yàn)結(jié)果(中子發(fā)射后)[51-53]進(jìn)行比較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 總結(jié)與展望

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

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