溫稠密物質(zhì)物態(tài)方程的理論研究

馬桂存張其黎 宋紅州 李瓊 朱希睿 孟續(xù)軍

(北京應(yīng)用物理與計算數(shù)學(xué)研究所,北京 100089)(2016年11月23日收到;2016年12月26日收到修改稿)

專題:高壓下物質(zhì)的新結(jié)構(gòu)與新性質(zhì)研究進(jìn)展

溫稠密物質(zhì)物態(tài)方程的理論研究

馬桂存?張其黎 宋紅州 李瓊 朱希睿 孟續(xù)軍

(北京應(yīng)用物理與計算數(shù)學(xué)研究所,北京 100089)(2016年11月23日收到;2016年12月26日收到修改稿)

本文詳細(xì)地介紹了溫稠密物質(zhì)物態(tài)方程的理論模型,其中包括液體變分微擾理論、化學(xué)圖像模型、離化電離平衡模型、平均原子模型和INFERNO模型;給出了混合物物態(tài)方程的計算方法;對第一原理分子動力學(xué)和量子蒙特卡羅方法進(jìn)行了介紹;對一些典型材料(如氫、氘、氦、氙、金、鎢等)在溫稠密區(qū)的物態(tài)方程進(jìn)行了計算和總結(jié);分析了離解、電離效應(yīng)對物態(tài)方程的影響.

溫稠密物質(zhì),物態(tài)方程,等離子體

1 引 言

材料的物態(tài)方程是反映物質(zhì)的壓強、密度、溫度之間的函數(shù)關(guān)系的方程[1].它在工程和科學(xué)技術(shù)中有著廣泛的應(yīng)用背景,如天體物理中的恒星構(gòu)造、地球物理中的地核運動以及核聚變等領(lǐng)域中都需要知道材料的物態(tài)方程.在恒星體的內(nèi)部(如太陽以及銀河系中的其他類星體中),其物質(zhì)的狀態(tài)實際上是一種高溫、高密度的等離子體狀態(tài).地球的內(nèi)核實際上也是一種熔融的溫稠密的等離子體態(tài).

隨著科學(xué)技術(shù)的快速發(fā)展,目前人們能夠在實驗室中實現(xiàn)材料的高溫高壓狀態(tài),如激光慣性約束聚變和磁約束聚變.美國的國家點火裝置、中國的神光裝置以及即將在歐洲建造的國際熱核聚變實驗堆裝置等都是實現(xiàn)材料在極高溫極高壓狀態(tài)的大型設(shè)備.在這些研究領(lǐng)域中所涉及的材料的物態(tài)方程其溫度密度范圍都非常寬廣,從室溫到幾千萬度,壓強從常壓到幾千萬大氣壓.在如此寬廣的溫度密度變化范圍內(nèi),物質(zhì)經(jīng)歷著固態(tài)、液態(tài)、氣態(tài)、等離子體態(tài)等各種聚集態(tài)的變化,發(fā)生著相變、離解、電離等各種復(fù)雜的物理過程.在寬廣的溫度與密度范圍內(nèi),如何準(zhǔn)確給出材料的狀態(tài)方程是一項富有挑戰(zhàn)性的科學(xué)課題.

在高溫高壓下,物質(zhì)全部熔化,出現(xiàn)部分離化或電離等復(fù)雜的物理變化過程.這種部分離化的強耦合電子離子系統(tǒng)給理論描述帶來了很大的困難.托馬斯-費米理論只適用于極高壓強的物質(zhì)狀態(tài)(通常為幾千萬甚至億萬大氣壓).該模型不考慮原子的殼層結(jié)構(gòu),認(rèn)為物質(zhì)中的原子處于完全電離的狀態(tài).在千萬大氣壓的壓強段,物質(zhì)中的原子實際上是處于部分離化的狀態(tài),即一部分電子仍然被束縛在原子核周圍,另一部分電子則處于游離狀態(tài).也就是說,物質(zhì)在高溫高壓“過渡區(qū)”的狀態(tài)實際上是一種部分離化的等離子體狀態(tài).由于在這個區(qū)域內(nèi),粒子之間的相互作用很強,多體微擾理論不適用.而且,由于存在多種粒子的激發(fā)態(tài),使得統(tǒng)計物理計算十分困難.

用于該區(qū)域的理論模型有物理圖像法和化學(xué)圖像法.物理圖像法中的逸度展開方法[2]只能用于高溫低密度系統(tǒng).對于溫稠密系統(tǒng)的物態(tài)方程,一般都采用化學(xué)圖像方法計算.化學(xué)圖像方法需要先給出系統(tǒng)的自由能,如俄羅斯科學(xué)家Fortov等[3]提出的自由能模型,以及美國和歐洲科學(xué)家Saumon,Chabrier,van Horn提出的自由能模型,即SCVH模型[4].液體微擾理論和粒子之間的相互作用勢是該類模型的理論基礎(chǔ)[5?7].在用SCVH模型對氫高壓物態(tài)方程的研究中,他們發(fā)現(xiàn)了一種等離子體相變,即氫從導(dǎo)電率很低的絕緣態(tài)到導(dǎo)電率很高的金屬態(tài)的轉(zhuǎn)變.

除了上述物理與化學(xué)圖像外,處理部分離化的等離子體系統(tǒng)還可以采用經(jīng)驗?zāi)Ｐ?如Kerley[8]提出的離化電離平衡模型.將該模型與液體微擾理論相結(jié)合,可以獲得寬廣的溫度與密度范圍內(nèi)的物態(tài)方程[9].該模型是建立在經(jīng)驗內(nèi)插基礎(chǔ)之上,它有很多的經(jīng)驗參數(shù),確定起來很困難.有些參數(shù)需要靠沖擊波實驗以及其他一些實驗才能確定[10].

平均原子模型是將等離子體中的復(fù)雜離子系統(tǒng)等效為一個原子,用一個原子的性質(zhì)來表達(dá)整個系統(tǒng)的性質(zhì). 最有代表性的就是Rozsnyai[11]的EOSTA模型和Liberman[12]的INFERNO模型.該模型在過渡區(qū)物態(tài)方程的計算中也取得了相當(dāng)?shù)某晒?朱希睿和孟續(xù)軍[13]對Rozsnyai的EOSTA模型進(jìn)行了推廣,并將其用于Au,Al,Hg等材料的電子熱能熱壓計算.通過改進(jìn)計算方法,INFERNO模型已經(jīng)用于Al,Be,U等材料高壓雨貢紐的計算[14,15].Barshalom和Oreg[16]通過修改維里定理的表達(dá)方式,得到了與實驗符合得較好的冷能冷壓.最近,馬桂存等[17]又將該模型用于金高壓雨貢紐的計算.上述平均原子模型的共同特點是考慮了束縛態(tài)電子的壓致離化和溫致離化效應(yīng),能夠比較好地用于過渡區(qū).但該模型也有很多不足之處,主要是沒有考慮離子本身的熱運動,以及沒有考慮其他離子的熱運動對所研究的平均原子內(nèi)的電子態(tài)的關(guān)聯(lián)影響.

在物態(tài)方程的數(shù)值模擬方面,目前主要采用量子分子動力學(xué)和量子蒙特卡羅方法.分子動力學(xué)(如VASP[18])只能用于較低的溫度,一般認(rèn)為其適用的溫度不超過10 eV.盡管人們也提出了高溫數(shù)值模擬方法,如無軌道分子動力學(xué)方法[19].但是,如何構(gòu)造密度泛函和粒子之間的相互作用勢仍然是沒有解決的問題.量子蒙特卡羅方法雖然是用于多體相互作用系統(tǒng)比較好的數(shù)值模擬方法.但由于其計算量大,計算方法復(fù)雜,目前僅用于低原子序數(shù)材料的物態(tài)方程計算[20,21].

混合物的物態(tài)方程在實際應(yīng)用中也非常重要,如人們在研究木星的構(gòu)造時,就需要氫和氦混合物的物態(tài)方程.無論是用物理圖像法,還是化學(xué)圖像法計算混合物的物態(tài)方程,其計算量將成倍地增加.此外,如何給出不同種類粒子之間的相互作用勢也是一個非常復(fù)雜的問題.目前的第一原理計算只能模擬上百個原子,所以在目前的計算條件下,完全用第一原理方法來給出混合物的物態(tài)方程也是不現(xiàn)實的.因此,尋找從單質(zhì)物態(tài)方程構(gòu)造混合物物態(tài)方程的經(jīng)驗方法,就成為構(gòu)造混合物物態(tài)方程的主要方法.其中,體積相加模型就是一種獲得混合物物態(tài)方程的有效方法[22].

本文從兩方面來介紹溫稠密物質(zhì)物態(tài)方程的計算方法.在第2節(jié)中,主要介紹溫稠密物質(zhì)物態(tài)方程的理論模型,其中,包括考慮原子內(nèi)部電子熱激發(fā)效應(yīng)的流體變分微擾理論、化學(xué)圖像模型SCVH、離化電離平衡模型IEQ、平均原子模型和INFERNO模型以及混合物態(tài)方程的計算方法.在第3節(jié)中,主要介紹第一原理分子動力學(xué)方法和量子蒙特卡羅方法,以及這兩種方法在氫、氘、氦、金、鎢等材料物態(tài)方程計算中的應(yīng)用結(jié)果.最后,給出本文的總結(jié).

2 物態(tài)方程的理論模型

2.1 液體變分微擾理論

液體是常見的一種物質(zhì)狀態(tài),是固體熔化或氣體在低溫下被壓縮形成的一種狀態(tài),它的物態(tài)方程在物態(tài)方程的理論研究中占有十分重要的地位.由于液體分子在空間排列無序,以及粒子間存在很強的相互作用,理論上給出液體狀態(tài)方程都采用變分微擾論方法.在給定粒子之間相互作用勢的基礎(chǔ)上,以硬球模型或軟球模型為參考態(tài)來建立計算液體物態(tài)方程的理論計算方法,被稱之為液體變分微擾論(FVT)[23].傳統(tǒng)的液體物態(tài)方程理論模型是不考慮系統(tǒng)中粒子內(nèi)部電子的激發(fā).這種近似在低溫下是可行的,但隨著溫度的升高,電子被激發(fā)到高量子態(tài),甚至電離.在本節(jié)中我們將考慮這種熱效應(yīng)對物態(tài)方程的貢獻(xiàn),并稱考慮熱效應(yīng)的液體變分微擾論為IFVT.

對于給定溫度和密度的熱力學(xué)系統(tǒng),亥姆霍茲自由能是系統(tǒng)的熱力學(xué)特征函數(shù),完全決定系統(tǒng)的熱力學(xué)性質(zhì).所以,模型的建立從構(gòu)造系統(tǒng)的自由能開始.相互作用系統(tǒng)的自由能可以分為三個部分:

其中,Fid是理想氣體部分,Fnon-id是相互作用引起的自由能,Fint是內(nèi)部運動的自由能.根據(jù)液體變分微擾,系統(tǒng)的自由能為

其中,參數(shù)ε=243.1 K,α=12.8,ra=4.42?,W=1.28?;近距離的非物理吸引勢用指數(shù)衰減排斥勢代替,A,B的取值保證函數(shù)值和一次導(dǎo)數(shù)在銜接處連續(xù).對上述自由能變分,可求得與自由能極小值相對應(yīng)的硬球直徑.

考慮到熱效應(yīng)修正,自由能需要增加原子的殼層電子躍遷項Fint,

其中,gα和εα分別是氙原子第α個能級的簡并度和能量,數(shù)值取自NIST數(shù)據(jù)庫[24].

基于上述液體變分微擾方法,我們計算了氙的雨貢紐壓強和溫度,并考查了熱效應(yīng)修正的影響,其結(jié)果如圖1和圖2所示,其中FVT和IFVT分別為傳統(tǒng)的液體變分微擾方法和考慮了熱效應(yīng)修正的液體變分微擾方法.從6 g/cm3處開始熱效應(yīng)修正顯著降低了雨貢紐壓強和溫度,使其與實驗數(shù)據(jù)[25,26]符合更好.這是由于殼層電子躍遷能夠吸收沖擊波能量,使得原子的熱運動能量減少,從而顯著降低溫度和壓強.由此可以看出,采用原子的殼層電子躍遷項能夠比較準(zhǔn)確地描述液體變分微擾方法的熱效應(yīng)修正,使得雨貢紐壓強和溫度的計算值在較大范圍內(nèi)得到顯著改善,與傳統(tǒng)的液體變分微擾方法相比,熱效應(yīng)修正使得模型的適用范圍得到了較大擴(kuò)展.

圖1 氙的Hugoniot壓強隨密度的變化,沖擊初態(tài)為165 K,2.96 g/cm3Fig.1.The Hugoniot of xenon,the initial state is 165 K,2.96 g/cm3.

圖2 氙的Hugoniot溫度隨密度的變化,沖擊初態(tài)同圖1Fig.2.The Hugoniot temperatu re changes with density,the initial state is the same as that in Fig.1.

2.2 化學(xué)平衡模型

當(dāng)溫度逐漸升高或密度進(jìn)一步增大,就會引起分子的離解和原子的電離,形成分子、原子、離子和電子的混合物.對于這樣一個復(fù)雜的系統(tǒng),計算該系統(tǒng)物態(tài)方程的方法主要有自由能極小方法.Saumon,Chabrier和van Horn提出了一種計算系統(tǒng)自由能的方法,即SCVH模型[4?6].該模型是一種化學(xué)圖像法,即事先假定系統(tǒng)中包含有哪幾種物質(zhì)成分,如分子、原子、離子以及電子等,然后根據(jù)系統(tǒng)的自由能極小來確定各種成分的比例.下面我們將以氫為例來具體介紹SCVH模型.

在低溫下,氫以分子形式存在,當(dāng)溫度升高或密度增加時,會發(fā)生氫分子的離解和氫原子的電離.SCVH模型分兩步給出了系統(tǒng)的自由能,首先給出氫分子和氫原子混合物的自由能,然后再給出氫離子(即質(zhì)子)和電子組成的等離子體系統(tǒng)的自由能.對于氫分子和氫原子組成的混合物,它的自由能為

其中,Fid是分子或原子的質(zhì)心運動動能;Fconf表示相互作用引起的自由能,即H2-H2,H-H以及H2-H之間相互作用導(dǎo)致的自由能;Fint是由分子內(nèi)部自由度和原子內(nèi)部自由度導(dǎo)致的自由能,即分子的振動與轉(zhuǎn)動以及原子內(nèi)電子熱運動的能量;Fqm是考慮分子和原子自身運動的量子效應(yīng)修正,它是一個小量,由W igner[27]和Kirkwood[28]的?2展開方法給出.對于相互作用部分Fconf,則采用以硬球系統(tǒng)為參考系的Weeks-Chand ler-Andersen形式的液體微擾論方法[29?31]計算.準(zhǔn)確給出分子與分子、原子與原子以及分子與原子之間的相互作用勢是本模型的關(guān)鍵.H2分子之間的有效兩體相互作用勢用沖擊波實驗數(shù)據(jù)擬合確定,H原子之間以及H2分子與H原子之間的相互作用勢由第一原理方法給出.由于SCVH模型涉及壓致離化區(qū)域,所以Fint部分采用Hummer和Mihalas[32]提出的占據(jù)概率方法計算,

式中,下標(biāo)i=1,2分別對應(yīng)于H2,H;ωiα表示束縛態(tài)α存在的概率(0≤ ωiα≤ 1);Ni表示系統(tǒng)中分子或原子的數(shù)目;εiα,giα分別表示分子或原子內(nèi)部運動的能級和簡并度.

對于由質(zhì)子和電子組成的等離子體系統(tǒng),SCVH模型采用屏蔽的一元等離子體模型描述,即認(rèn)為質(zhì)子在由電子組成的負(fù)電荷背景中運動,離子之間的相互作用勢用屏蔽的庫侖勢Veff描述.這種模型被稱之為SCOP(screened one component plasma),它的自由能為[6]

所以,整個系統(tǒng)的自由能為

在給定的溫度和密度下,通過求系統(tǒng)的自由能極小,可以給出各種粒子成分的含量.此外,通過熱力學(xué)關(guān)系式,可以求出系統(tǒng)的內(nèi)能、壓強等物理量.圖3給出了氫原子的濃度隨溫度密度的變化關(guān)系[5].溫度從2884—18197 K在對數(shù)坐標(biāo)下等間隔分布.從圖3可以看出,在低溫(T＜3000 K)下,氫分子含量比較多,氫原子占的比例比較少.隨著溫度的升高,氫原子的濃度在低密度區(qū)由小變大逐漸接近于1,這表明溫致離解效應(yīng)在該區(qū)域起作用.其次,隨著密度的增加氫原子濃度在減小,在密度為0.5—1.0 g/cm3內(nèi)有一個明顯的轉(zhuǎn)變,它反映的是壓強導(dǎo)致的離解和離化效應(yīng).即在給定的溫度下,當(dāng)密度大于1.0以后,系統(tǒng)中的氫分子將不復(fù)存在,取而代之的是氫原子或氫離子.圖4給出了氫的物態(tài)方程,即壓強隨溫度密度的變化[4].其中的兩條低溫等溫線出現(xiàn)了明顯的轉(zhuǎn)折,它反映了在該區(qū)域出現(xiàn)了等離子體相變(plasma phase transition,PPT),即從導(dǎo)電率很低的絕緣體狀態(tài)轉(zhuǎn)變到導(dǎo)電率很高的金屬狀態(tài).該相變類似于氣液相變,即物理量存在不連續(xù)變化,且有臨界點存在.這種相變的內(nèi)在機(jī)理就是壓致離化,即由壓強或密度的增加而導(dǎo)致的凝聚體物性的突變.關(guān)于PPT的詳細(xì)討論請看參考文獻(xiàn)[6].

圖3 在一系列的溫度下氫原子的濃度隨密度的變化Fig.3.The concentration of Hatomchanges with density at a series of temperatu res.

圖4 氫在一系列溫度下的等溫線Fig.4.The equation of state of hyd rogen at a series of temperatu res.

2.3 離化電離平衡模型

上節(jié)中的SCVH模型依賴于粒子之間的詳細(xì)相互作用,給出粒子之間的相互作用勢是一項復(fù)雜而繁重的工作.在本節(jié)中,我們將在經(jīng)驗?zāi)Ｐ偷幕A(chǔ)上,給出氘的物態(tài)方程.在給定的溫度和密度下,對于由單一成分組成的系統(tǒng),其能量和壓強通常寫成如下三項式形式:

其中,Ec,Pc為冷能、冷壓;Eion,Pion為離子的熱能、熱壓;Ee,Pe為電子的熱能、熱壓.對于由氘分子,或氘原子組成的系統(tǒng),其冷能、冷壓都用如下經(jīng)驗公式描述[9]：

式中,EB是結(jié)合能;δ= ρ/ρ0;α,ν是可調(diào)參數(shù),可根據(jù)實驗數(shù)據(jù)給出.

離子部分的貢獻(xiàn)由修正的硬球模型(即CRIS模型[34?36])給出:

式中N〈?〉0是一級勢能項[34];為硬球系統(tǒng)的堆積因子(packing fraction of hard sphere);?E1,?E2,?P1,?P2是高價修正項[35,36].

熱激發(fā)的電子對物態(tài)方程的貢獻(xiàn)采用離化電離平衡模型(IEQ)計算[8].設(shè)在體積為V和溫度為T的系統(tǒng)中,含有N個粒子,它包括不帶電的中性原子和帶正電的離子.原子的核電荷數(shù)為Z,原子量為W.每個原子或離子的電子結(jié)構(gòu)是通過指定每個軌道上的電子占據(jù)數(shù)來確定,然后,根據(jù)量子力學(xué)原理來確定每個電子的軌道能量,系統(tǒng)的宏觀量是各種離子成分貢獻(xiàn)的統(tǒng)計平均.在描述電子運動時,要區(qū)分電子的兩種狀態(tài),一種是束縛態(tài),一種是巡游態(tài).處于束縛態(tài)的電子只能圍繞原子核運動,一般用原子軌道波函數(shù)描述;處于巡游態(tài)的電子可以在整個系統(tǒng)中運動,一般用平面波描述.該模型的優(yōu)點是它考慮了各種離子成分對物態(tài)方程的貢獻(xiàn).系統(tǒng)的自由能可以表示為[8]

其中,qe為系統(tǒng)的配分函數(shù),qz為帶電荷為z的離子的配分函數(shù).(21)式中的求和是對各種離子態(tài)進(jìn)行的,各種離子態(tài)的電荷數(shù)分別為z=0,1,2,3,...,Z,共包含Z+1種離子態(tài);uz為帶電荷z的離子的基態(tài)能量;a(z)為一個自由電子的自由能,電中性要求每個離子平均含有z個自由電子;δz表示每個原子的電荷漲落.(22)式中的求和表示對每個離子的能級進(jìn)行求和,其中εz(n),gz(n)分別表示離子的各個量子態(tài)的能量和相應(yīng)的簡并度.為了能夠使上述模型準(zhǔn)確地描述電子的電離態(tài),必須做一些修正,如連續(xù)譜能量修正、電荷漲落效應(yīng)和體積漲落效應(yīng)修正等.

根據(jù)上述物態(tài)方程理論模型,我們計算了液氘在一定的溫度與密度范圍內(nèi)的物態(tài)方程,并計算了液氘在一定的壓力范圍內(nèi)的雨貢紐.圖5和圖6給出了液氘的低壓和高壓雨貢紐,同時也給出了雨貢紐實驗點、第一原理計算以及各種理論模型的結(jié)果.

圖5 液氘的雨貢紐(低壓部分)Fig.5.The Hugoniot of liquid deuterium(lowpressure part).

圖6 液氘的雨貢紐(高壓部分)Fig.6.The Hugoniot of liquid deuterium(high pressu re part).

目前用于測量液氘雨貢紐的實驗有四類:第一類是化爆和二級輕氣炮[37],它們給出的雨貢紐壓強只有20 GPa(見圖中Gas gun);第二類是球面內(nèi)爆[38],它的壓強達(dá)到110 GPa(見圖中Explosive)；第三類是Z-pinch[39],它達(dá)到的壓強也有100 GPa;第四類是強激光實驗[40?42],它達(dá)到的壓強有350 GPa(見圖中Laser).較早期的強激光實驗[40,41](見圖中Silva和Collins)由于其測量的方法導(dǎo)致的誤差非常大.所以,它給出的雨貢紐線最大壓縮比也非常大.后期的Z-pinch、球面內(nèi)爆以及強激光實驗結(jié)果[42]都給出了較小的最大壓縮比.最早的第一原理計算結(jié)果是由美國科學(xué)家Militzer和Ceperley[43]用量子蒙特卡羅方法在2000年給出(見圖中QMC).后來的第一原理計算結(jié)果,如采用直接積分的量子蒙特卡羅方法[44]以及第一原理分子動力學(xué)計算結(jié)果[45](見圖中QMD)都表明,液氘雨貢紐的最大壓縮比為4.65左右.同時,在這一時期,物態(tài)方程的理論模型也得到了很大的改進(jìn).Ross[46]的線性混合模型提出得比較早,該模型給出的雨貢紐曲線的壓縮比也比較大,與早期的激光實驗結(jié)果差不多.后來提出的物態(tài)方程理論模型,如液體變分微擾理論[7,23](見圖中FVT),Saumon和Chabrier[4?6]的化學(xué)平衡模型(見圖中SCVH),以及在本節(jié)中以Kerley的Sesame物態(tài)方程模型為基礎(chǔ)新做的物態(tài)方程(見圖中Newmodel[9]),都給出了較小的壓縮比.

2.4 INFERNO模型

INFERNO模型是在固體凝膠模型中加入一個原子而形成的平均原子模型.在通常的固體結(jié)構(gòu)中,帶正電的離子按照一定的方式排列,在每個離子周圍又有很多電子,它們代表核芯電子與價電子(或巡游電子).核芯電子只能在束縛它的原子核周圍運動,而巡游電子屬于整個系統(tǒng),可以在物質(zhì)內(nèi)部自由運動.在量子統(tǒng)計自洽場INFERNO模型中,認(rèn)為正離子是均勻分布在整個空間的,同時,為了保持系統(tǒng)的電中性,也有均勻分布的同樣多的負(fù)電荷存在,這種模型被稱為凝膠模型.現(xiàn)在,在凝膠模型中加入一個原子,該原子的核電荷位于系統(tǒng)中央,為了保持系統(tǒng)的電中性,在以原子核為中心,原子半徑為大小的空間內(nèi),所有的正電荷將被排斥掉,剩下的是按照量子力學(xué)規(guī)律分布的負(fù)電荷.這種模型通常被稱之為INFERNO模型[12,47].

在INFERNO模型中,球內(nèi)電子的波函數(shù)滿足相對論形式的單電子狄拉克方程[12,47]:

式中α,β是Dirac矩陣,c是光速.在原子球內(nèi)(或muffi n-tin球內(nèi))的勢函數(shù)V(r)可以表示為

其中,Z是核電荷數(shù),ρ(r)是球內(nèi)電子電荷密度分布,R是原子球半徑,υ為拉格朗日乘子.mu ffi n-tin球外的勢函數(shù)可以表示為

在INFERNO模型中,電子的狀態(tài)被分為束縛態(tài)、共振態(tài)、巡游態(tài)(或散射態(tài))三種類型,電子的波函數(shù)要通過自洽求解量子力學(xué)的Dirac方程獲得.電子所受到的勢場是一個自洽的多體勢,它包含了電荷之間的庫侖能,以及多體效應(yīng)引起的交換關(guān)聯(lián)能.共振態(tài)、散射態(tài)的出現(xiàn)給理論計算帶來了很大的困難[47],特別是在高能量、高角動量的情況下,電子波函數(shù)要發(fā)生劇烈振蕩,這些問題的出現(xiàn)使理論計算非常復(fù)雜.為了克服上述困難,美國利弗莫爾實驗室的科學(xué)家提出了用位相-振幅形式的波函數(shù)方法[14]來描述電子態(tài),克服了高能量的連續(xù)譜在計算上困難,得到了很好的結(jié)果.2009年,法國原子能科學(xué)研究中心的科學(xué)家提出了采用大規(guī)模并行計算的方法[15]來克服高能量積分,高角動量求和困難.2011年,俄羅斯科學(xué)家也用類似的方法計算了INFERNO模型的物理量[48].除了在計算方法上得到改進(jìn)外,最近,以色列科學(xué)家又對INFERNO模型在壓強計算上進(jìn)行了改進(jìn)[16,49,50],提出了根據(jù)維里定理計算壓強的方法,得到了與實驗一致的很好結(jié)果.

我們利用INFERNO模型,對Au在寬廣的溫度與密度范圍內(nèi)的物態(tài)方程進(jìn)行了計算[17].圖7給出了金屬Au中的電子壓強隨密度的變化,每一條線對應(yīng)一個溫度.溫度從log10(2.3×103)到log10(4.8×106),等間距分布.從圖7可以看出,在比較低的溫度下,電子系統(tǒng)的壓強隨密度都發(fā)生了非單調(diào)變化.引起這個變化的主要原因是Au原子中的外殼層電子的狀態(tài)發(fā)生了變化,即從束縛態(tài)變成了巡游態(tài),這種轉(zhuǎn)變被稱之為絕緣體金屬轉(zhuǎn)變或Mott轉(zhuǎn)變[51].對發(fā)生這種轉(zhuǎn)變的理解并不是非常困難.設(shè)想初始時,系統(tǒng)中的Au原子處于密度小于10 g/cm3的某個狀態(tài),顯然,在這種狀態(tài)下,Au的外殼層電子被原子核束縛,整個系統(tǒng)是絕緣體.如果在溫度不變的條件下,對系統(tǒng)進(jìn)行壓縮,在某個密度下,Au原子的外殼層電子波函數(shù)就要發(fā)生重疊,它將引起外殼層的d電子和s電子從束縛態(tài)變成巡游態(tài).轉(zhuǎn)變完成,整個系統(tǒng)將變成金屬.這一轉(zhuǎn)變的主要特征是其d電子能譜由寬變窄,再變寬,即d電子由束縛態(tài)變?yōu)楣舱駪B(tài),再變?yōu)檠灿螒B(tài).這種現(xiàn)象被稱之為壓致電離效應(yīng)[47,51],即由壓強引起的電子態(tài)的轉(zhuǎn)變.在用INFERNO模型對Au全區(qū)域物態(tài)方程的計算過程中,在某些溫度密度點處經(jīng)常會遇到共振能態(tài),它是由壓強或溫度變化引起的電子狀態(tài)的變化.共振能態(tài)的主要特點是其態(tài)密度有尖銳峰出現(xiàn),而且峰寬非常窄,這給計算帶來了很大的困難.為了克服這一困難,我們對共振能態(tài)的特點進(jìn)行了仔細(xì)的分析,得到了計算連續(xù)譜的波函數(shù)的有效方法[47],這種方法對處理共振能態(tài)特別有效.

對于離子熱運動部分,我們采用一種內(nèi)插模型來描述固態(tài)或液態(tài)系統(tǒng)的物態(tài)方程,它就是大家常用的自由體積模型[1].在給出了物態(tài)方程的各個部分貢獻(xiàn)之后,就可以按照三項式形式的物態(tài)方程計算材料的雨貢紐.圖8給出了Au的理論雨貢紐.為了比較,我們也在圖8中給出了利用一些加載設(shè)備如化爆、二級輕氣炮等從實驗上獲得的金材料的高壓雨貢紐[52?55].從圖8可以看出理論與實驗符合得非常好.在很高的壓強下材料的雨貢紐曲線發(fā)生了彎曲,它是由原子的內(nèi)殼層電子的離化從而吸收沖擊波能量所致.

圖7 在一系列溫度下Au中的電子熱壓隨密度的變化Fig.7.The electronic pressure of gold changes with density at a series of temperatures.

圖8 Au的高壓雨貢紐曲線,其中,實線是理論結(jié)果Fig.8.The Hugoniot curve of gold at high pressu re,solid line is the resu lt of theory.

2.5 平均原子模型

在前面幾節(jié),我們已經(jīng)涉及了物質(zhì)中的原子的離化或電離效應(yīng),如在SCVH模型、INFERNO模型以及IEQ模型中,都考慮到了原子的電離效應(yīng).在每一種模型中,都給出了用一定的方法來處理周圍的帶電粒子對所考慮的離子的影響.INFERNO模型對束縛電子與自由電子之間的轉(zhuǎn)換做了自洽處理.實際上,在部分離化的等離子體區(qū)域,還有另一種經(jīng)常使用的平均原子模型,即Rozsnyai[11]給出的平均原子模型.它是基于溫稠密物質(zhì)的等效原子模型思想,即認(rèn)為在整個系統(tǒng)中每個原子都是等價的.這樣每原子周圍的環(huán)境都是一樣的,因而導(dǎo)致在原子邊界上的電荷密度是連續(xù)的,周期性邊界條件要求電荷密度的導(dǎo)數(shù)也應(yīng)連續(xù).這些條件界定了離子內(nèi)部束縛電子的占住方式.自由電子和準(zhǔn)自由電子的存在,以及離子內(nèi)部電子細(xì)致組態(tài)的考慮,使模型在計算上變得比較復(fù)雜.在Hartree-Fock理論的框架下,電子的運動方程為

其中 Vsc(r)為自洽勢,?i(r)為單電子波函數(shù),εi為單電子能級.在球坐標(biāo)系下分離變量,可得到徑向波函數(shù)Pnl(r)滿足的方程.考慮到電荷密度及其導(dǎo)數(shù)的連續(xù)性,徑向波函數(shù)Pnl(r)應(yīng)滿足下列兩種形式的邊界條件：

第一種邊界條件取為

周期場中的電子,其波函數(shù)一般要受到簡諧波的調(diào)制.對應(yīng)來講,波函數(shù)按第一類邊界條件銜接時是因為整個物質(zhì)內(nèi)部處于最低能態(tài),按第二類邊界條件銜接意味著物質(zhì)內(nèi)部已處于最高能態(tài).由于物質(zhì)內(nèi)部調(diào)制波的周期是最基本周期的整數(shù)倍,那么原子間波函數(shù)的連接就不只是第一類和第二類邊界條件,還可能是中間過渡的邊界條件.

嚴(yán)格確定原子的波函數(shù)取第一類邊界條件還是第二類邊界條件是沒有太大意義的.當(dāng)原子中的電子波函數(shù)取第一類邊界條件與第二類邊界條件之間的允許值時,電子的能級就具有了一定的范圍,這個能級范圍也可以叫做“能帶”.不過這個能帶的起源與固體理論中能帶的起源本質(zhì)是不同的.這個“能帶”對確定有界原子束縛態(tài)波函數(shù)的能態(tài)密度具有重要意義.能帶的出現(xiàn)使得高溫稠密等離子體的原子結(jié)構(gòu)求解更加復(fù)雜.為了簡化模型,一般取第一類邊界條件.

原子胞中的電子被分為束縛電子和自由電子,其數(shù)密度可以寫為

自由電子的數(shù)密度ρf(r)本應(yīng)由自由電子的正能態(tài)波函數(shù)來確定,但由于正能態(tài)能量、軌道量子數(shù)都是無限的,為了快速計算,采用分波加統(tǒng)計的方法解決,并且具有較高的精度.自由電子數(shù)密度ρf(r)可以寫為

其中,V(r)是自洽勢,μ是化學(xué)勢,T是溫度,l0是分波的界限,其取值視所需計算量和精度所定.束縛態(tài)電子的數(shù)密度ρb(r)為

關(guān)于壓強的計算,由于已經(jīng)把一部分波函數(shù)能延伸到原子邊界處的束縛電子劃歸準(zhǔn)自由電子用分波法處理,所以剩下的束縛電子的徑向波函數(shù)就基本不能延伸到原子邊界處,這樣電子壓強就只能與非束縛電子有關(guān).這樣我們得到的電子壓強就可以只由三部分組成:

其中,Pk為自由電子產(chǎn)生的動壓,Pr為分波電子產(chǎn)生的壓強,Pex+corr為交換、關(guān)聯(lián)效應(yīng)綜合產(chǎn)生的負(fù)壓.

在含溫有界原子模型下,精確的原子總能量由四部分組成:束縛電子的組態(tài)平均能量Ebav,自由電子的單電子能Ef、束縛電子與自由電子間的相互作用能Eb-f,自由電子與自由電子間的相互作用能E f-f,即

關(guān)于壓強和能量的各部分貢獻(xiàn)的詳細(xì)計算請看參考文獻(xiàn)[13,56].由于我們的含溫有界模型是自動包含溫度、密度效應(yīng)的,并且由于部分分波法的引入大大降低了高溫下原子參數(shù)和狀態(tài)方程參量的計算量,因而我們可以在保證相當(dāng)?shù)木认驴焖俚玫饺我鉁囟让芏认碌臓顟B(tài)方程參數(shù).

圖9和圖10展示了本模型計算的Au在溫度分別為0.1,5.0,10.0,50,100,500,1000 eV下,電子壓強P和原子能量E隨物質(zhì)密度的大范圍變化(密度為1.0—5000 g/cm3,基態(tài)能E0=?17865.2288 a.u.).從計算結(jié)果中可以看出,在密度為20 g/cm3以及70 g/cm3處均顯示出了較強的原子殼層效應(yīng),這是由于原子殼層的壓制離化所導(dǎo)致的.在圖10中,原子的能量不再隨密度的增加而穩(wěn)步上升,甚至?xí)霈F(xiàn)隨密度的增加原子能量降低的現(xiàn)象.這是由于本模型考慮的是孤立原子的能量,在密度比較稀時,低殼層會出現(xiàn)一些空軌道,而密度的增加標(biāo)志著壓力的增大,會把一些外殼層的電子壓回低殼層的空軌道上,使得電子的總勢能下降,從而造成整個原子的能量會有所降低.

圖9 在一系列的溫度下Au電子壓強隨密度的變化Fig.9.The electronic pressu re of gold changes with density at a series of temperatu res.

圖10 在一系列的溫度下Au原子能量隨密度的變化Fig.10.The energy per atomof gold changes with density at a series of temperatu res.

2.6 混合物物態(tài)方程

在激光慣性約束聚變以及天體物理的行星構(gòu)造研究中,人們經(jīng)常會遇到混合物情況.如何給出混合物物態(tài)方程是備受關(guān)注的一個問題.例如在研究木星構(gòu)造時,就要知道氫與氦混合物的狀態(tài)方程.木星表層中[57]氫占74.2%,氦占23.1%,其他元素占0.027%.前面幾節(jié)介紹的物態(tài)方程計算方法都是對單一成分而言的,對于由多種元素組成的混合物,如果按照液體微擾論方法計算物態(tài)方程,那就需要給出不同粒子之間的相互作用勢,這在實際問題中很難做到.如果按照第一原理計算方法求物態(tài)方程,對于不同的混合比都一一去做數(shù)值計算,其計算量也是相當(dāng)大的.所以,有必要尋找一種從單質(zhì)物態(tài)方程計算混合物物態(tài)方程的經(jīng)驗方法.線性混合模型或體積相加模型可以比較好地給出混合物的物態(tài)方程.在這個模型中,認(rèn)為在等溫等壓條件下,混合物的體積是單質(zhì)成分在與混合物有相同的壓強和溫度下的體積之和,即

其中,N是混合物中組元的總數(shù),ni是第i組元的摩爾分?jǐn)?shù).如果用質(zhì)量百分比χi來表示,即

其中,ρ是混合物的密度,ρi是組元i的密度.系統(tǒng)的內(nèi)能為

系統(tǒng)的熵為

其中,Ei,Si分別是組元i的內(nèi)能和熵;Smix是混合熵增.

下面以氫與氦混合物為例,具體討論體積相加模型在物態(tài)方程計算中的應(yīng)用.圖11給出了用體積相加模型計算的氫氦混合物的物態(tài)方程,氫與氦的原子個數(shù)比為NH:NHe=2,同時我們將第一原理分子動力學(xué)(QMD)計算結(jié)果[22]也標(biāo)注在圖中,可以看出兩者符合得比較好,說明體積相加模型是計算混合物物態(tài)方程的一個很好的近似計算方法.但是,用第一性原理分子動力學(xué)和量子蒙特卡羅方法對氫氦混合物物態(tài)方程的詳細(xì)研究表明:在某些溫度密度范圍內(nèi),線性混合模型與第一原理結(jié)果相比卻給出了高達(dá)12%誤差[22].我們也用分子動力學(xué)方法計算了氘氦混合物的物態(tài)方程[58].定義混合比x=NHe/(NH+NHe).圖12給出了氘氦混合物在壓強P=10 GPa,不同溫度、不同混合比下的體積差:?Vmix/V=(VLM?V)/V隨混合比x的變化,其中V是實際體積,VLM是混合模型給出的體積.從圖12可以看出,在高溫下,混合比x接近0.5時,偏差最大,最高可達(dá)到8%.研究表明,氦混合到氫中會對氫分子的離解產(chǎn)生一定的影響.同時它也對氫原子的電離產(chǎn)生一定的影響.這種影響在某些情況下是無法用線性混合模型來表達(dá)的.

圖11 在不同溫度下氫氦混合物的物態(tài)方程,其中,實線是體積相加模型的結(jié)果,數(shù)據(jù)點是第一原理分子動力學(xué)結(jié)果[22]Fig.11.The equation of state for Hand He mixing at d iff erent temperatu res.lines are the resu lts of add itive volume,dots fromQMD[22].

圖12 在壓強為10 GPa時不同溫度下體積差隨混合比的變化Fig.12.The volume d iff erence changes with mixing fraction of Heliumat pressu re 10 GPa.

如果我們定義系統(tǒng)的過剩壓強為總壓強與理想氣體壓強之差,即Pex=P?Pid,則在線性混合的假設(shè)下,讓每個單一成分按照等溫等Pex混合,其混合體積仍然是各成分的體積之和,那么按照這種方式得到的混合物物態(tài)方程與第一原理的差別會更小[59].從熱力學(xué)上看,這種讓過剩壓強相等的混合模型相當(dāng)于認(rèn)為混合系統(tǒng)是由獨立的平均有效原子組成,各種有效原子之間不存在相互作用,混合系統(tǒng)的過剩自由能就是每個有效原子的自由能之和.

3 物態(tài)方程的第一原理計算

在溫稠密物質(zhì)系統(tǒng)內(nèi)部,粒子之間的關(guān)聯(lián)相互作用非常強,隨著溫度和密度的變化在溫稠密物質(zhì)內(nèi)部會發(fā)生部分離化、強關(guān)聯(lián)以及量子效應(yīng)等一系列復(fù)雜過程,特別是在金屬絕緣體轉(zhuǎn)變區(qū)域,很難找到一種有效的理論方法來描述其變化行為.第一原理分子動力學(xué)和量子蒙特卡羅方法為我們提供了認(rèn)識該區(qū)域物態(tài)變化行為的有效數(shù)值模擬工具,兩者都是解決多體問題的有效方法.第一原理分子動力學(xué)是建立在密度泛函理論基礎(chǔ)上的有限溫度數(shù)值模擬方法,量子蒙特卡羅方法是建立在路徑積分基礎(chǔ)上的一種隨機(jī)行走方法.

3.1 第一原理分子動力學(xué)方法

計算材料在有限溫度下狀態(tài)方程的數(shù)值方法,目前比較常用的是基于密度泛函理論的分子動力學(xué)方法,通常稱為量子分子動力學(xué)(QMD)方法[18,60],該方法通過Born-Oppenhaimer近似,將電子離子問題部分解耦合,每個離子都在電子與其他離子形成的有效勢場中運動,離子受到的力通過Hellman-Feynman定理給出,離子運動按照經(jīng)典力學(xué)規(guī)律描述.在離子運動的每個時刻中,電子運動由求解該時刻離子位置構(gòu)型下的電子Hamiltonian系統(tǒng)來確定.QMD方法嚴(yán)重依賴贗勢和交換關(guān)聯(lián)勢的準(zhǔn)確性,同時隨著溫度升高,需要計算的能帶數(shù)目增加,導(dǎo)致計算量大大增加,因此QMD不適合在很高溫度下的材料模擬.在本節(jié)中,我們將介紹利用QMD方法對氦、金、鎢的物態(tài)方程以及鎢的電學(xué)性質(zhì)的計算.

在分子動力學(xué)模擬中,一般采用超原胞和周期性邊界條件,模擬的原子數(shù)越多,所需的計算量就越大,通常用64或128個原子就可以達(dá)到理想的結(jié)果.電子和離子之間的相互作用勢選用超軟贗勢或投影綴加波贗勢,交換關(guān)聯(lián)泛函選用廣義梯度近似Perdew-Burke-Ernzerhof[61,62]形式,倒空間布里淵區(qū)積分采用Monkhorst-Pack[63]方法或用單Gamma點方法;平面波展開的波函數(shù)截斷能量一般為幾百eV.平面波截斷能量和k點的數(shù)目在實際計算時都要經(jīng)過優(yōu)化,以達(dá)到滿足實際精度的要求.離子系統(tǒng)的運動釆用NVT系綜描述,由Nose-Hoover[64]恒溫器來控制離子溫度,電子溫度分布則根據(jù)Fermi-Dirac統(tǒng)計得到.

3.1.1 氦的物態(tài)方程

我們用分子動力學(xué)VASP程序計算了氦在10000 K和30000 K的等溫線(見文獻(xiàn)[58]圖1).計算表明,在溫度不超過30000 K,密度不超過5 g/cm3的范圍內(nèi),用化學(xué)模型來計算狀態(tài)方程是合理的.在更高的溫度密度區(qū)域,則需要進(jìn)一步改進(jìn)模型.圖13顯示了根據(jù)物態(tài)方程計算得到的氦的一次沖擊和二次沖擊Hugoniot曲線,同時還顯示了Nellis等[65]的實驗數(shù)據(jù).從圖13可以看出,QMD和化學(xué)平衡模型(CM)的計算結(jié)果符合得很好,且一次沖擊Hugoniot與實驗數(shù)據(jù)一致,二次沖擊Hugoniot結(jié)果與實驗數(shù)據(jù)相比稍微偏軟,這與實驗點的精度有關(guān).盡管如此,理論仍然在實驗的誤差范圍.

圖13 氦的一次和二次沖擊Hugoniot曲線Fig.13.The Hugoniot of heliumfromfi rst and second shock.

3.1.2 金的物態(tài)方程

我們用分子動力學(xué)方法計算了Au密度在19.3—42 g/cm3,溫度在1000—50000 K區(qū)間的物態(tài)方程,并在此基礎(chǔ)上得到了Hugoniot曲線[66],如圖14所示[52,67?70].圖中D是沖擊波速度,u是粒子速度,其中的數(shù)據(jù)點是實驗結(jié)果,粉色的曲線是馬桂存等[67]在2005年用QEOS方法計算的結(jié)果.在低壓區(qū)(u＜1.5 km/s),我們的結(jié)果與QEOS基本一致,且與早期的實驗結(jié)果符合得很好.在1.5 km/s＜u＜6 km/s區(qū)間,我們的結(jié)果與QEOS方法計算的結(jié)果差別較大,Yokoo研究組[52]的實驗結(jié)果與我們的計算結(jié)果更接近,QEOS方法計算的結(jié)果偏軟,其原因是QEOS方法在做參數(shù)擬合時只用到了早期的化爆和氣炮實驗結(jié)果.在更高壓區(qū),目前的實驗手段只有激光加載才能達(dá)到,但由于實驗中各種復(fù)雜的因素導(dǎo)致激光加載實驗數(shù)據(jù)的精度較差、數(shù)據(jù)點分散,因此實驗結(jié)果不能用來標(biāo)定理論數(shù)據(jù).由此可以得出結(jié)論:僅利用低壓的實驗數(shù)據(jù)擬合參數(shù)得到的QEOS理論狀態(tài)方程,不能保證在高壓時的準(zhǔn)確性;為了能使其用于慣性約束聚變數(shù)值模擬中,需要更多的高精度的高壓實驗數(shù)據(jù)來標(biāo)定理論物態(tài)方程的參數(shù).

圖14 金的沖擊波速度與粒子速度關(guān)系Fig.14.The shock velocity changes with particle’s velocity for gold.

3.1.3 鎢的物態(tài)方程

從理論上說,數(shù)值模擬膨脹態(tài)金屬流體系統(tǒng)的金屬-非金屬轉(zhuǎn)變非常困難,由于系統(tǒng)在極端條件下的無序性和熱激發(fā)行為,體系成分復(fù)雜,包括離子、自由電子、中性原子、分子及團(tuán)簇等,且各成分高度瞬時變化;此外物理過程也非常復(fù)雜,包括各種分解、激發(fā)、電離及其逆過程,體系成分及其中發(fā)生的物理過程強烈地受周圍環(huán)境的影響,必須考慮多體效應(yīng)以及關(guān)聯(lián)相互作用.基于有限溫度密度泛函理論的量子分子動力學(xué)方法提供了有效模擬極端情況下流體特性的工具.該方法對離子實采用經(jīng)典分子動力學(xué)描述,對電子采用有限溫度量子處理,并結(jié)合線性響應(yīng)理論的Kubo-Greenwood公式,在統(tǒng)一的框架下計算材料的物態(tài)方程、光學(xué)性質(zhì)、電學(xué)性質(zhì)和熱力學(xué)性質(zhì).

我們采用第一原理分子動力學(xué)方法研究了膨脹鎢在較寬的溫度與密度范圍的物態(tài)方程,并與已有的實驗結(jié)果進(jìn)行了比較,如圖15所示.以溫度為300 K時的體積V0作為參考體積.其中,用QMD計算的等溫線的溫度分別為0,300,500,1000,1500 K.在圖15中,我們還給出了Konstantin等[71]采用WC壓砧,MgO-Au作為壓標(biāo),給出的一系列溫度下的壓強實驗數(shù)據(jù).理論與實驗符合得很好.

圖15 金屬鎢在膨脹區(qū)的物態(tài)方程Fig.15.The equation of state for tungsten in the expanded region.

3.1.4 鎢的電導(dǎo)率

我們利用分子動力學(xué)VASP程序計算體系的電導(dǎo)率,主要通過Kubo-Greenwood公式[72,73]計算電導(dǎo)的實部隨頻率的變化:

式中Ψi,k和Ψj,k為Kohn-Sham方程本征態(tài),εi,k和εj,k為對應(yīng)的本征值;f(εi,k)是電子占據(jù)數(shù);ω(k)是k點權(quán)重.

圖16 在不同的溫度下鎢的靜態(tài)電導(dǎo)率隨密度變化行為,其中,在密度為3 g/cm3附近,電導(dǎo)率出現(xiàn)了反轉(zhuǎn),表明出現(xiàn)金屬-非金屬轉(zhuǎn)變Fig.16. The static conductivity of W changes with density at d iff erent temperatures. Near density 3 g/cm3,the behavior of conductivity changed abruptly,which means that metal tonon-metal transition occu rred.

金屬與非金屬區(qū)別在于在電子能譜中是否有能帶間隙,而帶隙的出現(xiàn)可以通過電導(dǎo)反映出來.由于離化引起熱激發(fā),體系直流電導(dǎo)將隨溫度的增加而增大,但是隨著密度的減小又將導(dǎo)致電導(dǎo)率的降低.體系的實際電導(dǎo)隨溫度和體積的變化關(guān)系是上述兩種行為的競爭結(jié)果,體系發(fā)生金屬-非金屬轉(zhuǎn)變需要低電子濃度以及達(dá)到一定程度的無序性,最終導(dǎo)致電子的局域化并引起金屬-非金屬相變.

3.2 量子蒙特卡羅方法

路徑積分蒙特卡羅(PIMC)方法是比較適合的模擬系統(tǒng)有限溫度性質(zhì)的第一性原理方法[76].該方法已經(jīng)很好地用于描述Bose系統(tǒng)的熱力學(xué)性質(zhì),但對于Fermi系統(tǒng),由于波函數(shù)的交換反對稱性導(dǎo)致所謂的“負(fù)符號”問題限制了它的應(yīng)用.為了解決這個問題,Ceperley等提出了約束路徑積分蒙特卡羅(RPIMC)方法[76?78],Hu等[21]利用RPIMC方法計算了一套氘的狀態(tài)方程數(shù)據(jù).Filinov等[79?81]認(rèn)為RPIMC方法引入了不可控的近似,他們發(fā)展了直接費米路徑積分蒙特卡羅(DPIMC)方法,并成功計算了非理想氫等離子體平衡態(tài)的熱力學(xué)性質(zhì).我們采用DPIMC方法對氫的狀態(tài)方程進(jìn)行了計算,并與QMD和RPIMC的結(jié)果進(jìn)行了比較.

3.2.1 DPIMC方法簡介

眾所周知,量子系統(tǒng)的熱力學(xué)性質(zhì)完全由其配分函數(shù)Z決定.對于由Ne個電子和Np個質(zhì)子組成的二元混合系統(tǒng),其配分函數(shù)Z可表示為

其中,β=1/kBT;q=q1,q2,...,qNp和r={r1,r2,...,rNe}分別表示質(zhì)子和電子的坐標(biāo);σ={σ1,σ2,...,σNe}表示電子的自旋態(tài).ρ是密度矩陣,可以表示為路徑積分的標(biāo)準(zhǔn)形式:

根據(jù)統(tǒng)計物理學(xué)的熱力學(xué)公式,可以得到系統(tǒng)能量和壓強的表達(dá)式.能量的計算公式為βE=?β?ln Q/?β,壓強公式為βp=?ln Q/?V=[α/3V?lnQ/?α]α=1.以上計算狀態(tài)方程的公式中都包含多重積分,不能直接進(jìn)行計算,但很適合采用標(biāo)準(zhǔn)的蒙特卡羅方法進(jìn)行數(shù)值計算.

3.2.2 DPIMC計算結(jié)果

我們對氫等離子體狀態(tài)方程進(jìn)行了大量的計算.圖17是125000 K溫度下的壓強隨密度的變化,我們將計算結(jié)果與其他數(shù)值方法進(jìn)行了比較,圖中的QMD表示基于密度泛函理論的量子分子動力學(xué)方法[82],RPIMC表示約束路徑積分蒙特卡羅方法[77].在低密度下,關(guān)聯(lián)效應(yīng)和簡并效應(yīng)都較弱,三種數(shù)值方法計算的結(jié)果符合得很好.但在高密度下,關(guān)聯(lián)效應(yīng)和簡并效應(yīng)都很強,三種方法計算的結(jié)果存在比較大的偏差.QMD計算的壓強比兩種路徑積分方法更高,可能的原因有以下兩點:1)在高溫下,QMD一般采用無軌道形式,其電子動能算符的密度泛函形式采用Thomas-Fermi模型,該模型高估電子的壓強;2)在密度泛函理論中,交換-關(guān)聯(lián)勢通常采用對零溫電子氣的蒙特卡羅模擬結(jié)果進(jìn)行擬合得到的經(jīng)驗形式,其能否準(zhǔn)確描述高溫電子的交換-關(guān)聯(lián)效應(yīng),沒有經(jīng)過檢驗.兩種路徑積分的結(jié)果也存在偏差,與它們在處理“負(fù)符號”問題時所引入的近似有關(guān).

圖17 溫度T =125000 K時壓強與密度的關(guān)系(1 bar=105Pa)Fig.17.The pressu re changeswith density at temperature T=125000 K(1 bar=105Pa).

由于QMD與DPIMC計算能量的方式不同,它們計算的能量不能直接比較,在圖18中我們比較了DPIMC和RPIMC計算的能量.在大部分溫度范圍內(nèi),兩者都符合得很好,只在低溫(T=31250 K)下有明顯差別,如圖18的內(nèi)插圖中最左邊的點所示.在該點,Γ=3.37,Θ=0.48,關(guān)聯(lián)和簡并效應(yīng)都不能忽略.這說明DPIMC和RPIMC對關(guān)聯(lián)和簡并效應(yīng)考慮的近似程度是不一樣的,在目前的實驗技術(shù)條件下,在該溫度密度下有可能通過精密的實驗來進(jìn)行驗證,并為更準(zhǔn)確的理論模型提供參考.

圖18 電子數(shù)密度為0.05962 A?3時,能量與溫度的關(guān)系Fig.18.The energy changes with temperatu re at electron density 0.05962 A?3.

4 結(jié) 論

針對溫稠密物質(zhì)狀態(tài)的變化特點,我們介紹了一系列能夠計算該區(qū)域物態(tài)方程的物理模型.建立在液體微擾論基礎(chǔ)上的化學(xué)圖像模型SCVH能夠比較好地描述物質(zhì)從分子到原子、從原子到離子的離解和電離過程.該模型的不足之處是它依賴于粒子之間的詳細(xì)相互作用勢.建立在經(jīng)驗?zāi)Ｐ突A(chǔ)之上的離化電離平衡模型(IEQ),也是一種近似的化學(xué)平衡模型.它也能夠比較好地給出物質(zhì)在過渡區(qū)的物態(tài)方程.其缺點是該模型有一些經(jīng)驗參數(shù)需要由實驗來確定.平均原子模型在概念上雖然簡單,但它卻是描述過渡區(qū)原子電離的很好模型.在INFERNO模型中,由于對共振電子態(tài)的特殊處理,從而避免了電子電離帶來的物理量跳躍.這是INFERNO模型的最大優(yōu)點.但是該模型也有不足之處,主要是它沒有考慮離子的熱運動,以及沒有考慮其他離子對該離子的關(guān)聯(lián)影響.

第一原理分子動力學(xué)和量子蒙特卡羅方法是計算溫稠密區(qū)物態(tài)方程的有效方法.這兩種數(shù)值模擬方法的最大特點是它能夠在量子力學(xué)的框架內(nèi)高效、準(zhǔn)確地處理分子離解、原子電離等復(fù)雜的物理過程,而且,還能夠?qū)爻砻芪镔|(zhì)的電學(xué)性質(zhì)、光學(xué)性質(zhì)等輸運特性做很好的計算.兩者共同特點是模擬的粒子數(shù)少,計算量大.量子分子動力學(xué)方法向高溫推廣也有一定的困難.量子蒙特卡羅方法向多電子束縛態(tài)的低溫推廣也有一定的困難.

等離子體相變(即PPT)是在溫稠密區(qū)發(fā)現(xiàn)的一個有趣的物理現(xiàn)象,但是在一些分子動力學(xué)模擬計算中卻沒有發(fā)現(xiàn)這一現(xiàn)象[83,84],所以有必要從理論和實驗上對該問題做進(jìn)一步地探討.其次,在溫稠密區(qū)出現(xiàn)的多元素混合物的相分離現(xiàn)象(如氫氦混合物[85])也是值得進(jìn)一步研究的課題.這些現(xiàn)象將直接影響人們對行星(如木星)構(gòu)造的認(rèn)識.

[1]Xu X S,Zhang W X 1986In troduction tothe Theory of Equation of State(Beijing:Science Press)(in Chinese)[徐錫申,張萬箱 1986實用物態(tài)方程理論導(dǎo)引(北京:科學(xué)出版社)]

[2]Rogers F J,Young D A1997Phys.Rev.E56 5876

[3]Fortov V,Iakubov I,Kh rapak A2006Physics of Strongly Coupled P lasma(Oxford:Ox ford University Press)

[4]Saumon D,Chabrier G,van Horn HM1995Astrophys.J.Suppl.Ser.99 713

[5]Saumon D,Chabrier G 1991Phys.Rev.A44 5122

[6]Saumon D,Chabrier G 1992Phys.Rev.A46 2084

[7]Ju ranek H,Redmer R 2002J.Chem.Phys.117 1768

[8]Kerley G I1986J.Chem.Phys.85 5228

[9]Kerley G I2003Equations of State for Hydrogen and DeuteriumSandia National Laboratories Technical Report No.SAND 2003-3613

[10]Jin F Q 1999In troduction tothe Experimen ts of the Equation of State(Beijing:Science Press)(in Chinese)[經(jīng)福謙1999實驗物態(tài)方程導(dǎo)引(北京:科學(xué)出版社)]

[11]Rozsnyai BF 1972Phys.Rev.A5 1137

[12]Liberman D A1979Phys.Rev.B20 4981

[13]Zhu X R,Meng X J 2011Acta Phys.Sin.60 093103(in Chinese)[朱希睿,孟續(xù)軍 2011物理學(xué)報 60 093103]

[14]W ilson B,Sonnad V,Sterne P,IsaacsW 2006J.Quan t.Spectrosc.Radiat.Transfer99 658

[15]Penicaud M2009J.Phys.:Condens.Matter21 095409

[16]BarshalomA,Oreg J 2009High Energy Density Physics5 196

[17]Ma G C,Zhang Q L,Lu G 2017Chin.J.High Press.Phys.31 1(in Chinese)[馬桂存,張其黎,盧果 2017高壓物理學(xué)報31 1]

[18]Kress G,Fu rthmu ller J 1996Phys.Rev.B54 11169

[19]Lambert F,C lerouin J,Zerah G 2006Phys.Rev.E73 016403

[20]D river KP,Militzer B2012Phys.Rev.Lett.108 115502

[21]Hu S X,Militzer B,Goncharov V N,Skupsky S 2011Phys.Rev.B84 224109

[22]Vorberger J,Tamb lyn I,Militzer B,Bonev S A2007Phys.Rev.B75 024206

[23]Ju ranek H,Redmer R 2000J.Chem.Phys.112 3780

[24]Kramida A,RalchenkoYu,Reader J,N ISTASD TeamN ISTAtomic Spectra Database(ver.5.3)http://physics.nist.gov/asd[2015-1-2]

[25]NellisW J,van Thiel M,Mitchell AC 1982Phys.Rev.Lett.48 816

[26]Urlin V D,Mochalov MA,Mikhailova OL 1992High Press.Res.8 595

[27]W igner E 1932Phys.Rev.40 749

[28]Kirkwoord J 1933Phys.Rev.44 31

[29]Andersen HC,Chand ler D 1970J.Chem.Phys.53 547

[30]W eeks J D,Chand ler D,Andersen HC 1971J.Chem.Phys.54 5237

[31]W eeks J D,Chand ler D,Andersen HC 1971J.Chem.Phys.55 5422

[32]Hummer D G,Mihalas D 1988Astrophys.J.331 794

[33]Ebeling W,Forster A,Richert W,Hess H1988Physica A150 159

[34]Kerley G I1980J.Chem.Phys.73 469

[35]Kerley G I1980J.Chem.Phys.73 478

[36]Kerley G I1980J.Chem.Phys.73 487

[37]Nellis W J,Mitchell AC,van Thiel M,Devine G J,Trainor R J 1983J.Chem.Phys.79 1480

[38]Boriskov G V,Bykov AI,II’kaev R I,Selemir V D,Simakov G V,Trunin R F,Urlin V D,Shuikin AN,NellisW J 2005Phys.Rev.B71 092104

[39]Knudson MD,Hanson D L,Bailey J E,Hall C A,Asay J R 2001Phys.Rev.Lett.87 225501

[40]Da Silva L B,Celliers P,Collins GW,Budil KS,Holmes N C,Barbee TW,Hammel BA,Kilkenny J D,W allace R J,Ross M,Caub le R,Ng A,Chiu G L B1997Phys.Rev.Lett.78 483

[41]Collins G W,Da Silva L B,Celliers P,Gold D,Foord M,Wallace R,Ng A,W eber S,Budil K,Caub le R 1998Science281 1178

[42]Hicks D G,Boehly TR 2009Phys.Rev.B79 014112

[43]Militzer B,Ceperley D M2000Phys.Rev.Lett.85 1890

[44]Bezkrovniy V,Filinov V S,KrempD,Bonitz M,Sch langes M,Kraeft W D,Levashov P R,Fortov V E 2004Phys.Rev.E70 057401

[45]Caillabet L,Mazevet S,Loubey re P 2011Phys.Rev.B83 094101

[46]Ross M1998Phys.Rev.B58 669

[47]Ma G C,Qi J,W ang M2015Chin.J.Comput.Phys.32 361(in Chinese)[馬桂存,齊進(jìn),王敏 2015計算物理32 361]

[48]Novikov V G,Ovechkin AA2011Math.Models Comput.Simu lations3 290

[49]BarshalomA,Oreg J 2007High Energy Density Phys.3 12

[50]BarshalomA,Oreg J 2006J.Quan t.Spectrosc.Radiat.Transfer99 35

[51]Kerley G I1987In t.J.Impact.Eng.5 441

[52]YokooM,KawaiN,Nakamura KG,KondoK2008Appl.Phys.Lett.92 051901

[53]Marsh S P 1980Los Alamos Shock Hugoniot Data(Berkeley:University of California Press)

[54]Al’tshu ler L V,Bakanova AA,Dudoladov IP,Dynin E A,Trunin R F,Chekin BS 1981J.Appl.Mech.Tech.Phys.22 145

[55]Jones AH,IsbellW M,Maiden C J 1966J.Appl.Phys.37 3493

[56]Zhu X R,Meng X J,Tian MF 2008Acta Phys.Sin.57 4049(in Chinese)[朱希睿,孟續(xù)軍,田明峰2008物理學(xué)報57 4049]

[57]Militzer B2005J.Low.Temp.Phys.139 739

[58]Zhang Q L,Zhang G M,ZhaoY H,Liu HF 2015Acta Phys.Sin.64 094702(in Chinese)[張其黎,張弓木,趙艷紅,劉海風(fēng)2015物理學(xué)報64 094702]

[59]Danel J F,Kazand jian L 2015Phys.Rev.E91 013103

[60]Car R,ParrinelloM1985Phys.Rev.Lett.55 2471

[61]PerdewJP,Burke K,ErnzerhofM1996Phys.Rev.Lett.77 3865

[62]W ang Y,PerdewJ P 1991Phys.Rev.B44 13298

[63]Monkhorst HJ,Pack J D 1976Phys.Rev.B13 5188

[64]Nose S C 1984J.Chem.Phys.81 511

[65]NellisW J,Holmes N C,Mitchell AC,Trainor R J,GovernoG K,Ross M,Young D A1984Phys.Rev.Lett.53 1248

[66]Zhang Q L,ZhaoY H,Ma G C 2014Chin.J.High Press.Phys.28 18(in Chinese)[張其黎,趙艷紅,馬桂存2014高壓物理學(xué)報28 18]

[67]Ma G C,PeiW B,W ang M2005Annual Report of ICF272(in Chinese)[馬桂存,裴文兵,王敏 2005 ICF科技年報272]

[68]Koenig M,Faral B,Boudenne J,Batani D,Benuzzi A,Bossi S,Remond C,Perrine J P,Temporal M,Atzeni S 1995Phys.Rev.Lett.74 2260

[69]McQueen R,Marsh S 1960J.Appl.Phys.31 1253

[70]Al’tshu ler L,Krupnikov K,Brazhnik MI1958Sov.Phys.JETP7 614

[71]Konstantin D L,Pavel N G,Dorogokupets P I,Igor SS,Shatskiy A,Fei Y W,PashchenkoS V,Seryotkin Y V,HigoY J,Funakoshi K,Ohtani E 2013J.Appl.Phys.113 133505

[72]KuboR 1957J.Phys.Soc.Jpn.12 570

[73]G reenwood D A1958Proc.Phys.Soc.London715 585

[74]Mott N 1984Rep.Prog.Phys.47 909

[75]KorobenkoV N,Rakhel AD 2013Phys.Rev.B88 134203

[76]Ceperley D M1995Rev.Mod.Phys.67 279

[77]Pierleoni C,Ceperley D M,Bernu B,MagroW R 1994Phys.Rev.Lett.73 2145

[78]MagroW R,Ceperley D M,Pierleoni C,Bernu B1996Phys.Rev.Lett.76 1240

[79]Filinov V S,Bonitz M,Ebeling W,Fortov V E 2001P lasmas Phys.Con trolled Fusion43 743

[80]Filinov V S,Fortov V E,Bonitz M,KrempD 2000Phys.Lett.A274 228

[81]Filinov V S,Bonitz M,KrempD,Kraeft W D,Ebeling W,Levashov P R,Fortov V E 2001Con trib.P lasmas Phys.41 135

[82]W ang C,Zhang P 2013Phys.P lasmas20 092703

[83]Holst B,Redmer R,Desjarlais MP 2008Phys.Rev.B77 184201

[84]Redmer R,Ropke G 2010Contrib.Plasma Phys.50 970

[85]Militzer B2013Phys.Rev.B87 014202

PACS:64.10.+h,62.50.–p,91.60.Fe,52.25.KnDOI:10.7498/aps.66.036401

Theoretical study of the equation of state for warmdense matter

Ma Gui-Cun?Zhang Qi-Li Song Hong-Zhou LiQiong Zhu Xi-Rui Meng Xu-Jun

(Institu te of Applied Physics and Compu tational Mathematics,Beijing 100089,China)(Received 23 November 2016;revised manuscript received 26 December 2016)

In this paper,we present in detail various theoretical models for studying the equation of state of warmdense matter,including the fluid variational theory,the chemicalmodel,the ionization equilibriummodel,the average atommodel and INFERNOmodel.The method of calculating the equation of state of a mixture is alsogiven.The results fromthe fi rst principles molecu lar dynamics simu lation and the quantumMonte Carlosimu lation are alsoprovided.Typicalmaterials such as hydrogen,deuterium,helium,xenon,gold,tungsten,etc.are studied in warmdense region by using all themethods,showing the eff ects of dissociation and ionization in the equation of state.

warmdensematter,equation of state,plasma

10.7498/aps.66.036401

?通信作者.E-mail:ma_guicun@iapcm.ac.cn

?Corresponding author.E-mail:ma_guicun@iapcm.ac.cn