金陽 張平 李永軍 侯永? 曾交龍 袁建民2)
1) (國防科技大學文理學院物理系, 長沙 410073)
2) (中國工程物理研究院研究生院, 北京 100088)
在天體物理和慣性約束聚變研究中涉及到的溫稠密物質通常包含多種元素的混合, 并且每種元素還被電離成多種離子價態, 不同價態離子結構及其豐度將直接影響溫稠密物質的診斷及其物理性質.同時, 從電子結構計算出發來研究宏觀物理性質時, 還需要考慮溫度、密度效應對離子結構的影響.本文從不同價態離子的電子結構計算出發, 采用考慮了離子間相互作用的Saha 方程獲得了稠密環境下的離子豐度, 并使用超網鏈(hypernetted-chain)近似對鋁、金以及碳-氫混合物的徑向分布函數進行了計算, 結合離子周圍電子的密度分布, 最后獲得X-射線湯姆遜散射的彈性散射譜.在X-射線散射譜計算中, 計算了溫稠密物質中同時存在不同離子價態時的電子結構和徑向分布函數, 發現在相同的等離子體環境下不同價態離子的徑向分布函數和電子結構差別較大.這將對依賴于微觀統計過程的物理性質, 比如散射光譜, 將產生較大的影響.
過去幾十年, 由于激光技術的發展, 強激光不僅能夠在實驗室中產生溫稠密物質, 而且還可以產生X-射線為診斷溫稠密物質結構提供了手段, 使得高能量密度物理的研究得到快速發展.X-射線彈性散射譜能夠反映出溫稠密物質中的離子結構, 而離子結構的信息在計算中直接取決于離子的徑向分布函數和電子分布.同時, 徑向分布函數也能夠反映出等離子體中離子間的關聯效應[1?3].另外,通過離子結構的計算結果還可以直接獲得宏觀物理性質, 比如, 物態方程和離子輸運性質[4,5].所以,離子結構的精確計算對天體物理和慣性約束聚變[1]中所涉及到基本物理參數的研究有著重要應用價值.
在波恩-奧本海默近似下, 量子分子動力學(quantum molecular dynamics, QMD)[6?8]是在含溫密度泛函理論框架下計算電子結構, 通過能量變分獲得離子間的受力, 進而通過求解離子運動的牛頓方程可以很好地描述溫稠密物質中的電子、離子結構, 并能夠得到合理的宏觀物理性質.但在較高溫度、密度狀態下, 由于模擬系統中的高能量電子的描述需要更寬能量狀態區域, 這將使得計算量大大增加, 限制了高溫狀態下模擬系統中包含更多離子的可能性.為了快速的獲得溫稠密狀態下電子密度的空間分布, 人們使用Thomas-Fermi(TF)近似[9]直接計算離子周圍電子密度分布, 發展了無軌道的分子動力學(orbital free molecular dynamics,OFMD)[10]來計算高溫稠密物質的性質, 但對于束縛軌道的計算將帶來較大的誤差.在溫稠密狀態下, 由于溫度和密度的效應, 物質將會被電離成不同價態的離子, QMD 和OFMD 都不能有效的提取不同離子價態分布信息.對于傳統低密度等離子體中離子結構的計算, 通常忽略離子間的相互作用, 采用自由離子模型來計算離子的結構.但對于溫稠密物質, 離子間的耦合較強, 忽略稠密環境效應將會對離子結構計算產生較大影響.本課題組已經發展了多離子的分子動力學模型[11], 考慮了稠密環境下離子間相互作用對離子豐度影響, 采用分子動力學模擬給出離子的輸運性質.但在較高溫度下, 由于離子豐度分布較寬, 需要同時模擬更多的離子或者模擬更久的時間才能給出合理的計算結果, 大大增加了計算量.在已知粒子間相互作用勢的情況下, 為了能夠快速準確的獲得離子間的徑向分布函數, 通過耦合Ornstein-Zernike 方程, 超網鏈(hypernetted-chain, HNC)近似[12?14]可以給出合理的描述.
本文將首先使用計算原子結構的相對論程序FAC(flexible atomic code)[15]獲得不同價態離子的電子結構; 然后采用考慮了離子間相互作用的Saha 方程給出離子價態分布[14,16], 使用超網鏈近似來計算不同價離子間的徑向分布; 最后結合不同價離子的徑向分布函數以及離子周圍的電子分布,給出X-射線彈性散射譜隨角度變化關系.在此基礎上, 通過比較溫稠密區域不同價離子的徑向分布函數和X-射線的彈性散射譜, 本文討論了隨著等離子體溫度密度變化, 離子豐度對徑向分布函數的影響, 并與平均原子計算結果進行了對比, 研究不同統計過程對離子結構以及X-射線彈性散射譜的影響.
徑向分布函數在確定氣體、等離子體、流體等多粒子系統結構特性的研究中扮演了重要角色, 它是用統計方法來描述系統中粒子間相對位置變化.在類流體的描述中, 徑向分布函數和其他的關聯函數只依賴粒子之間的相對位置(r=|r1?r2|), 對于多種粒子組成的系統中兩個粒子之間的徑向分布函數可以表示為

其中,ri,rj表示粒子的位置,〈〉表示正則系綜的平均.na,nb和Na,Nb分別表示a,b類粒子數密度和粒子數.總關聯函數可以用徑向分布函數表示為hab=gab ?1.通過對總關聯函數的傅里葉變換,還可以獲得多組分混合物質的靜態結構因子:

δab在a,b為同類離子時為δab=1 ,a,b為不同種類離子時為δab=0.
在實際的數值計算中, 已知粒子間的相互作用時, 可以通過不同的方式獲得粒子的徑向分布函數, 比如分子動力學[10]、蒙特卡洛方法[17]以及HNC近似[12?14].為了快速獲得混合等離子體中的離子間徑向分布函數, 本文采用HNC 近似來計算溫稠密物質中離子間的關聯函數, 如下式:

其中,hab,cab分別為總關聯函數和直接關聯函數,Vab為a粒子和b粒子之間的勢能.在關聯函數的計算過程中, 還需要Ornstein-Zernike 方程形成閉合關系, 從而進行自洽迭代來求解.Ornstein-Zernike 方程表示如下:

其中nc表示c類粒子的離子數密度.
采用HNC 近似來計算粒子間的徑向分布函數, 首先需要獲得粒子間的相互作用勢.從頭計算粒子間的相互作用勢需要求解多體問題, 是非常困難的.人們通常建立不同的計算模型來描述離子間的相互作用, 比如Yukawa 模型、Deutsch 模型[18]、Klimontovich 和Kraeft(KK)模型[19]以及Kelbg模型[20]等.這些模型都是針對不同條件下提出計算離子間相互作用的模型, 都有一定的適用范圍.在溫稠密物質中, 考慮到存在不同價離子和自由電子, 以及電子、離子間的庫倫屏蔽效應, 在計算離子間的徑向分布函數時采用Yukawa 勢來描述離子間的相互作用:

其中kD為屏蔽系數, 是屏蔽長度的倒數.考慮到在溫稠密物質中離子間的相互作用勢主要被自由電子屏蔽, 所以, 在計算中屏蔽系數時只考慮自由電子對屏蔽的貢獻, 由公式:

在溫稠密物質的X-射線散射中, 將一束已知光譜強度的X-射線照射在溫稠密物質上, 在不同角度探測被散射的X-射線的光譜就可以獲取溫稠密物質的溫度、密度以及平均電離度的信息.散射X-射線的轉移動量依賴于散射角度:

其中,k表示在不同散射角上收集到X-射線動量和初始動量k0的差.散射譜的測量和物質的動態結構因子S(k,ω) 直接相關, 其中ω表示X-射線轉移到電子上的能量.根據Chihara 公式[21], 可以將動態結構因子近似分成彈性散射、非彈性散射部分,其中非彈性散射又可以簡單分為伴隨動量反沖或動能交換的自由電子散射和伴隨激發或者電離的束縛電子的散射.彈性散射部分是原子核以及隨原子核一起運動的電子對X-射線的散射, 對這部分散射可以忽略X-射線的能量損失, 彈性散射可以給出離子的結構信息.由于溫稠密物質中存在不同價態的離子, 所以在計算X-射線散射譜時需要考慮多種離子價態的貢獻.動態結構因子中彈性散射部分可以表示為形狀因子與靜態結構因子的乘積[22,23], 也被稱為Rayleigh 峰, 其在動量空間可以表達為

其中靜態結構因子由多組分的超網鏈近似獲得; 形狀因子由每價離子周圍分布的自由電子和束縛電子貢獻, 在傅里葉空間可以將形狀因子表示為

其中ρb(r) 為束縛電子密度, 根據組態的占據數乘上波函數模方求和得到:

其中Pn,Qn是求解Dirac 方程得到波函數的大、小分量,wn表示束縛軌道上的占據電子數目.自由電子密度由Thomas-Fermi 近似得到如下式所示:

其中p0(r)=[2V(r)c2+V2(r)]/c,μ為化學 勢, 由電中性條件得到,c為光速.

圖1 在溫度為104 K、離子數密度為1024 cm–3 時, Au1+和Au2+混合離子價態的徑向分布函數.Present 表示本文工作的計算結果(實線); HNC 表示文獻[2]中采用HNC 近似計算的結果(虛線 + 上三角); OCP 表示平均成一種價態離子的徑向分布函數; average 表示各價離子的徑向分布函數的算數平均Fig.1.The radial distribution function of Au1+, Au2+ in gold plasma with ion number density of 1024 cm–3 and temperature of 104 K.Present represents the result of the present paper (solid line); HNC represents the result of HNC approximation in the Ref.[2] (dotted line + upper triangle); OCP is the radial distribution function of mean charge-state ion; average labels the average results of the radial distribution function of different charge-state ions.
在溫稠密物質中由于多種價態離子同時存在,并且每一種價態離子所形成的勢場不同, 從而影響該離子周圍的其他離子分布, 徑向分布函數可以對離子的分布給出很好的描述.首先在溫度為104K、離子數密度為n= 1024cm–3的條件下, 使用HNC近似計算了金元素二價和一價等粒子數混合情況下的徑向分布函數, 如圖1 所示, 不同價離子間的相互作用勢都采用Yukawa 對勢形式進行計算.在圖中, 將計算結果與Wünsch 等[22]已發表的結果[2]進行了比較.因為在計算中都使用了Yukawa 勢函數, 并且都是使用HNC 近似方法計算得到的計算結果, 所以, 可以看出本文的計算結果和Wünsch等[22]結果基本上是重合的.從而也說明了, 發展的HNC 程序在計算不同價離子的關聯函數時是沒有問題的.為了進一步比較不同統計過程對徑向分布函數的影響, 首先按照離子豐度將不同價離子平均成一種價態(A= 1.5), 然后按照單組分形式可以給出一種平均離子的徑向分布函數(圖中OCP標記的計算結果); 也采用了多組分HNC 計算得到的不同價離子的徑向分布函數, 然后再根據豐度進行平均(圖中average 標記的計算結果).可以看出, 雖然都采用Yukawa 形式的勢函數, 但兩種統計的先后順序差異對計算結果產生了差別.這也說明了從微觀電子結構計算出發來獲得溫稠密物質的宏觀性質時, 不同的統計過程將影響所計算的宏觀物理性質.
為了進一步測試HNC 計算程序對于混合物質的計算, 在當前的理論框架下, 又計算了在溫度為2.0 × 104K、密度為2.5 × 1023cm–3CH 混合物質的徑向分布函數, 其中C4+和H+離子數密度比為1∶1.將計算結果和Wünsch 等[22]已發表的結果[2]進行比較, 如圖2 所示.從圖中可以看出, 計算的結果和Wünsch 等[22]結果基本一樣.從而進一步說明, 本文的HNC 程序能夠實現不同元素、同種元素不同價離子的徑向分布函數的計算.

圖2 在離子數密度為nH = nC = 2.5 × 1023 cm–3、溫度為T = 2 × 104 K 時, 碳-氫混合等離子體中C4+和H+的徑向分布函數; Present 是本文計算結果(實線), HNC 是文獻[2]中用HNC 近似的計算結果(虛線 + 上三角)Fig.2.The radial distribution function of C4+ and H+ in the CH mixture plasma at number density, nH = nC =2.5 × 1023 cm–3, and the temperature, T = 2 × 104 K.And solid lines label the present results; and dashed lines with upper triangles are the results of Ref.[2].
在計算離子間的相互作用勢時, 按照方程(5)的方式考慮了自由電子屏蔽效應, 構造了離子間的Yukawa 勢, 然后使用考慮了離子關聯效應的Saha 方程[10]計算了密度為8.1 g/cm3、不同溫度(10, 20, 40 和100 eV)下溫稠密鋁等離子體中離子豐度, 如表1 所示.從表中可以看出, 由于壓致電離效應, 在溫度為10 eV 時等離子體中只有3 種離子價態, 并且2 價、4 價離子的豐度很小, 主要是3 價離子, 這和常態下金屬鋁有3 個自由電子也相吻合.隨著溫度的升高, 2p 殼層被打開, 同時更多的電子電離導致離子周圍的勢阱變深, 不同價離子的能量差別越來越小, 所以, 即使在3 倍的固體密度下鋁等離子體中也同時存在更多價態的離子.
獲得了離子豐度和相互作用勢后, 就可以使用多組分的HNC 計算程序計算不同價態離子間的徑向分布函數.為了比較不同統計過程的計算結果, 也采用了平均原子模型+超網鏈近似的方法(AAHNC)[3]計算了平均離子的徑向分布函數.AAHNC 模型是首先將等離子體中的不同價態離子平均成一種價態, 在離子球里自洽計算電子的結構時考慮離子間的關聯效應, 其中的關聯函數也是由HNC 近似給出, 離子間的相互作用勢由修正的Gordon-Kim 方法計算.在計算不同價態離子的徑向分布函數時, 按照表1 給出的離子豐度, 在圖3中給出了鋁在密度為8.1 g/cm3、溫度分別為10,20, 40 和100 eV 時離子徑向分布函數的計算結果.圖中黃色實線是AAHNC 模型的計算結果, 不同顏色的虛線依次表示由多組分的HNC 近似得到的不同價態離子的徑向分布函數.從圖3 中可以看出, AAHNC 計算得到徑向分布函數都是位于各種不同價態離子計算結果的中間, 但由于平均電離度并不等于某個整數價態離子, 比如20、40 和100 eV的計算結果中, 它和中間價態離子的計算結果還是有差別的.另一方面, 即使在10 eV 時, AAHNC模型給出的平均電離度也是3 價, 但和直接采用多組分HNC 的計算的徑向分布函數也是有差別的,這就說明從不同的統計方法將會對離子結構的計算結果帶來影響.隨著溫度的升高, 平均電離度越來越大, 由于更多自由電子屏蔽, 相鄰價態離子的徑向分布函數差別越來越小.但由于離子價態分布范圍越來越寬, 平均價態離子計算結果和最高價、最低價離子計算結果的差別也越來越大.在這樣物理條件下, 一些非常依賴于微觀離子結構的宏觀物理性質, 比如光譜以及光譜的散射性質, 當從微觀離子結構計算出發, 使用不同統計過程將對計算結果產生較大影響.

表1 在密度為8.1 g/cm3、不同溫度下溫稠密鋁的離子豐度, A 表示平均電離度Table 1.Charge-state fractions of Al at density, 8.1 g/cm3, A is average charge state.
從電子結構計算出發來研究溫稠密物質體中的離子結構, 除了依賴于離子的空間分布, 還依賴于離子周圍的電子密度分布.所以, 在獲得了不同價離子豐度及徑向分布函數的基礎上, 還需要按照式(9)計算形狀因子(form factor).在溫稠密物質環境中, 由于離子間的強耦合以及電子簡并效應使得離子的高激發態較少, 因此在使用FAC 程序計算電子時, 只考慮了不同價態離子的基態電子結構[14].與圖3 一樣, 在圖4 中比較了不同價態離子和AAHNC 模型計算得到的形狀因子.由于不同價離子的中心對稱勢函數不同, 周圍電子密度分布的差異導致了不同價態離子的形狀因子差別, 特別對低k區域, 也就是外層電子的分布在小角度散射情況下對離子價態的依賴非常明顯.與徑向分布函數的計算結果類似, 雖然AAHNC 模型計算的形狀因子都是從不同價態離子的計算結果中間穿過,但隨著溫度的升高, 溫稠密物質中出現越來越多的不同價態離子, 特別對于100 eV 的計算結果, 由于很寬的離子價態分布導致最高價態、最低價態的計算結果差別很大.

圖3 在密度為8.1 g/cm3、不同溫度下溫稠密鋁中, 采用多組分HNC 近似給出不同價態離子(虛線)的徑向分布與AAHNC 模型(實線)計算結果的比較Fig.3.Different ion species pair distribution functions (dashed lines) of Al at density, 8.1 g/cm3, and different temperatures calculated by HNC approximation, comparing with that of AAHNC model (orange solid lines).

圖4 密度為8.1 g/cm3、不同溫度下溫稠密鋁等離子體中, 不同價態離子(虛線)、AAHNC 模型計算(實線)的形狀結構因子隨著散射角度的變化.Fig.4.Form factor of different ion species (dashed lines) of Al at density, 8.1 g/cm3, and different temperatures, comparing with that of AAHNC (orange solid lines).

圖5 在密度為8.1 g/cm3、溫度10 eV 下溫稠密鋁等離子體中, 離子結構隨不同散射角的變化關系.Multi-ion 表示本文使用的方法給出的計算結果(虛線); AAHNC 表示AAHNC模型的計算結果(實線)[3]; QMD 表示量子分子動力學計算結果(點-虛線)[24]; HNC 表示文獻[25]中采用HNC 近似計算離子結構的結果(點線); 帶有誤差范圍的點表示實驗結果[26]Fig.5.Ion feature for Al as function of k at a temperature of 10 eV and a density of 8.1 g/cm3: Multi-ion (dashed line), AAHNC (solid line)[3], QMD (dot-dashed line)[24],HNC (dot line)[25] and experimental data (points with error bars)[26].
為了更好地和實驗以及其他理論計算結果進行比較, 在獲得不同價離子的徑向分布和形狀因子基礎上, 通過方程(8)計算了溫度為10 eV、密度為8.1 g/cm3溫稠密鋁等離子體X-射線彈性散射隨不同散射角變化, 如圖5 所示.圖中的Multi-ion是表示本文使用的方法給出的計算結果.AAHNC是表示由一種平均價態離子的結構直接計算得到的結果.圖中的HNC 和AAHNC 方法類似, 都是采用超網鏈近似獲得離子間的關聯函數, 只是在計算電子-離子、離子-離子間的相互作用勢時采用了無規相近似(random-phase approximation, RPA)的方法.QMD 表示的量子分子動力學模擬的計算結果.從圖上可以看出, 除了在低k區域外, 考慮了多價態離子的計算結果和其他理論模型吻合都比較好.在低k區域的差別, 主要是由于在計算過程只考慮了離子對勢而忽略了離子的長程相互作用.對于稠密物質中, 由于自由電子的屏蔽, 忽略長程庫倫相互作用也是合理的近似.因此可以看出, 本文發展的多價態離子結果計算方法可以對溫稠密物質給出合理的描述.除了最高的兩個點外,理論計算結果和實驗吻合的都很好.在已發表的工作[27,28]中已經討論過, 如果認為該實驗測量對應的物質狀態處于非平衡狀態, 在理論模擬中考慮到電子、離子具有不同的溫度后, 可以得到了和實驗較吻合的結果.
本文從原子結構計算出發, 使用FAC 程序求解Dirac 方程獲得不同價態離子的電子結構及其能量; 根據自由電子分布, 采用Yukawa 模型構造了溫稠密物質中離子間相互作用勢函數; 在此基礎上, 通過考慮了離子間相互作用的Saha 方程計算了離子的豐度分布, 進而采用混合組分的超網鏈(HNC)近似計算溫稠密Au, Al 和CH 混合物質中不同價態離子徑向分布函數.計算結果與平均原子的計算結果進行比較, 發現統計的先后順序差異對離子結構及其電子分布的影響.基于不同價態的徑向分布函數及其周圍的電子分布, 按照Chihara 公式, 又計算了溫稠密Al 等離子體的 X-射線彈性散射譜, 雖然計算結果與其他理論和實驗結果吻合較好, 但不同價態離子的徑向分布函數和形狀因子差別較大, 進一步說明在溫稠密物質中考慮不同價態離子結構對宏觀物理性質的影響.