一個快速模擬熱稠密非平衡等離子體的碰撞輻射模型*

韓小英 李凌霄 戴振生 鄭無敵 古培俊 吳澤清

(北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所， 北京 100094)

在實(shí)驗(yàn)室和天體等離子體研究中， 原子/離子的激發(fā)、退激發(fā)以及電離、復(fù)合過程對等離子體的電離和能量平衡有著重要的影響.在激光等離子體作用的輻射流體模擬中， 需要在線計(jì)算等離子體的平均離化度和吸收/發(fā)射系數(shù).在現(xiàn)有的計(jì)算能力下， 通常采用比較簡單的平均原子(average atom， AA)模型進(jìn)行在線計(jì)算.隨著實(shí)驗(yàn)技術(shù)和計(jì)算能力的發(fā)展， 急需發(fā)展能夠在線計(jì)算的細(xì)致非平衡原子模型.本文報(bào)道了最新發(fā)展的多離化度平均離子碰撞輻射模型(multi-average ion collisional-radiative model， MAICRM).該模型用一個平均離子模擬等離子體中某一離化度所有離子的平均軌道占據(jù)數(shù)和布居， 即每個平均離子的軌道占據(jù)數(shù)為該離化度所有離子的軌道占據(jù)數(shù)的平均; 平均離子的布居等于該離化度離子的布居和.平均離子的軌道占據(jù)數(shù)和布居通過迭代求解速率方法得到.用該模型計(jì)算了Fe， Xe和Au非平衡等離子體的離化度分布， 計(jì)算結(jié)果與細(xì)致組態(tài)和超組態(tài)模型以及實(shí)驗(yàn)測量符合， 而計(jì)算量相對于細(xì)致組態(tài)/超組態(tài)大大降低.預(yù)期該方法能與輻射流體程序耦合， 實(shí)現(xiàn)細(xì)致非平衡原子模型的在線計(jì)算.

1 引 言

在實(shí)驗(yàn)室和天體等離子體的輻射流體模擬中，等離子體的離化度分布是其中關(guān)鍵的物理量[1，2].等離子體的電離平衡和離化度分布主要由原子(離子)的激發(fā)/退激發(fā)和電離/復(fù)合過程決定.當(dāng)密度高時(shí)， 等離子體接近局域熱動平衡狀態(tài)(local thermodynamic equilibrium， LTE)， 離子的軌道占據(jù)數(shù)和布居滿足Saha-Boltzmann分布; 當(dāng)密度低時(shí)， 等離子體處于冕區(qū)平衡狀態(tài)(coronal equilibrium， CE); 對于中等密度的等離子體， 則需要通過求解速率方程得到各離子的非平衡分布， 通常稱為碰撞-輻射(collisional-radiative， CR)模型.

目前模擬非平衡等離子體的方法主要有以下幾種(從考慮的能態(tài)由粗到細(xì)): 平均原子模型(average atom， AA)[3]、超組態(tài)方法(super-configuration accounting， SCA)[4]、細(xì)致組態(tài)方法(detailed configuration accounting， DCA)[5]、細(xì)致能級方法(detailed level accounting， DLA)[6].其中AA模型用一個平均原子模擬等離子體中各離子軌道的平均占據(jù)數(shù)給出平均離化度.AA模型的計(jì)算時(shí)間最短， 但是無法給出細(xì)致的離化度分布.DCA，SCA和DLA方法能給出大量能態(tài)(幾千或者幾萬個能態(tài))的細(xì)致分布， 但是計(jì)算時(shí)間比AA模型高幾個量級.此外， 上述細(xì)致原子模型對于不同溫度密度下的等離子體， 通常只能考慮有限的能態(tài)， 如何保證所選能態(tài)的完備性一直是該領(lǐng)域的一個難題.

在激光等離子體作用中， 等離子體要經(jīng)歷從低溫高密到高溫低密的寬范圍狀態(tài)， 目前通常采用AA模型在線計(jì)算等離子體的平均離化度和吸收、發(fā)射系數(shù).從吸收和發(fā)射總量來看， AA模型與實(shí)驗(yàn)測量符合[7，8].近年來， 隨著實(shí)驗(yàn)技術(shù)的發(fā)展， 人們不僅能測量等離子體的發(fā)射總量， 而且能夠進(jìn)行譜分辨的細(xì)致測量; 理論上也發(fā)展了所謂的“DCA”在線非平衡原子模型方法， 但是這里的“DCA”是用主量子數(shù)n刻畫組態(tài)[9]， 類似于n層次的超組態(tài)方法.文獻(xiàn)[10]比較了不同版“DCA”模型計(jì)算的發(fā)射總量， 發(fā)現(xiàn)計(jì)算結(jié)果有倍數(shù)差別.

針對細(xì)致、快速在線非平衡原子模型的需求，本研究組參考AA模型的物理圖像， 基于細(xì)致組態(tài)的速率方程， 發(fā)展了一種新的非平衡原子碰撞輻射模型， 即多離化度平均離子碰撞輻射模型(multiaverage ion collisional-radiative model， MAICRM).首先， 用一個平均離子來模擬某一離化度所有組態(tài)的軌道平均占據(jù)數(shù)， 該占據(jù)數(shù)通過求解軌道之間的激發(fā)和退激發(fā)速率方程得到.然后， 基于平均離子的軌道占據(jù)數(shù)， 通過求解離子之間的電離和復(fù)合速率方程， 得到各平均離子的布居.有了平均離子的軌道占據(jù)數(shù)和布居， 就可以得到等離子體的平均離化度和離化度分布以及細(xì)致的吸收和發(fā)射譜.該方法采用了兩個近似: 第一個近似認(rèn)為同一離化度離子的單軌道速率系數(shù)(不包含軌道的占據(jù)數(shù)和空穴數(shù))相等; 第二個近似將束縛軌道之間的躍遷(激發(fā)和退激發(fā)過程)和束縛軌道與連續(xù)軌道之間的躍遷(電離和復(fù)合過程)分兩步計(jì)算.

在DCA/SCA模型中， 每個離化度需要考慮幾百甚至幾千個組態(tài)/超組態(tài); AA模型只需考慮一個離化度的十個或者幾十個軌道之間的躍遷;MAICRM模型中， 最多只需考慮Z+1 個離子之間的電離、復(fù)合過程， 每個離子內(nèi)計(jì)算nmax個軌道之間的激發(fā)、退激發(fā)過程， 因此MAICRM相對于DCA/SCA模型能極大地減少計(jì)算量， 其計(jì)算量約為AA模型的Z+1 倍.

本文報(bào)道了用該方法計(jì)算的Fe， Xe和Au非平衡等離子體的離化度分布和平均離化度， 計(jì)算結(jié)果與DCA/SCA計(jì)算結(jié)果以及實(shí)驗(yàn)測量符合， 因此該方法有望在未來與輻射流體耦合， 實(shí)現(xiàn)在線細(xì)致非平衡原子模型計(jì)算.

2 理論方法

眾所周知， AA模型是用一個平均原子模擬所有離子軌道的平均占據(jù)， 該模型計(jì)算的平均離化度(核電荷數(shù)-束縛軌道占據(jù)數(shù))、總發(fā)射量等與實(shí)驗(yàn)測量基本相符， 但是無法獲得更細(xì)致的離化度分布等信息.從AA模型圖像出發(fā)， 考慮用一個平均離子模擬某個離化度所有離子的軌道平均占據(jù)數(shù)，從而可以得到更細(xì)致的離化度分布和不同離化度軌道的平均占據(jù)數(shù).基于該物理圖像， 可以從細(xì)致組態(tài)的速率方程出發(fā)， 通過對同一離化度所有組態(tài)的軌道占據(jù)數(shù)求平均， 得到各離化度離子的軌道平均占據(jù)數(shù)和布居的速率方程.

2.1 細(xì)致組態(tài)的速率方程

有ne個束縛電子的組態(tài)記為表示該組態(tài)各軌道的電子占據(jù)數(shù)， 組態(tài)的布居為PKne;gi和εi分別表示i軌道的統(tǒng)計(jì)權(quán)重和軌道能.該組態(tài)布居的速率方程如下:

式中， RE/D為激發(fā)/退激發(fā)速率系數(shù)， 包括電子碰撞和光過程.RI包括電子碰撞電離和光電離速率系數(shù)， RR包括輻射復(fù)合和三體復(fù)合速率系數(shù)，RAI為自電離速率系數(shù)， REC為電子俘獲速率系數(shù).各速率系數(shù)的計(jì)算公式詳見附錄A.注意， 這里的速率系數(shù) R 為單軌道速率系數(shù)， 不包含涉及軌道的占據(jù)數(shù)和空穴數(shù).

2.2 平均離子的速率方程

將(1)式兩邊乘以Kne組態(tài)的n軌道的占據(jù)數(shù)并對該離化度所有組態(tài)求和， 得到

求 和 過 程 中， 對 于u→l的 躍 遷， 利 用 了 關(guān) 系 式這里認(rèn)為同一離化度各組態(tài)的單電子軌道速率系數(shù) R 相同， 即 R 只與離化度有關(guān).該近似與SCA方法中認(rèn)為同一超組態(tài)的單軌道速率系數(shù)相同類似.

進(jìn)一步， 定義

其中，PΛne是所有具有ne個束縛電子的組態(tài)(Kne)的布居之和;是同一離化度所有組態(tài)(Kne)n軌道上占據(jù)數(shù)的平均值.據(jù)此， 定義一個平均離子它由以下三個物理量刻畫: 1)有ne個束縛電子; 2)布 居 為PΛne; 3)n軌 道 上 占 據(jù) 數(shù) 為為Λne離 子 的n軌 道 上 的 空 穴 數(shù).為Λne離子u和l軌道上占據(jù)數(shù)之積.將(3)式代入(2)式可得:

(4)式即為平均離子Λne的布居PΛne和軌道占據(jù)數(shù)的速率方程.在(4)式中是耦合的，很難直接通過(4)式求解考慮到一般況下， 電子的束縛-束縛躍遷速率比束縛-連續(xù)躍遷速率快， 軌道上電子占據(jù)數(shù)通常由激發(fā)和退激發(fā)過程決定， 而離化度的布居通常由電離和復(fù)合過程決定.因此， 將(4)式拆為下面兩個式子:(5)式為離子軌道占據(jù)數(shù)的速率方程， 包括軌道之間的激發(fā)和退激發(fā)過程.(6)式為離子布居的速率方程， 包括電離和復(fù)合過程.通過迭代求解(5)式和(6)式， 可以得到一組收斂的具體的計(jì)算過程為: 第一步， 根據(jù)給定的等離子體參數(shù)， 包括電子溫度(Te)、輻射場溫度(Tr)、離子溫度(Ti)和電子密度(Ne)， 計(jì)算激發(fā)/退激發(fā)速率系數(shù) RE/D; 第二步， 將激發(fā)/退激發(fā)速率系數(shù)代入(5)式， 假設(shè)每個離子初始占據(jù)數(shù)為基態(tài)占據(jù)， 再用隱式牛頓迭代法得到每個離子收斂的第三步， 根據(jù)每個離子的軌道占據(jù)數(shù)計(jì)算電離/復(fù)合速率系數(shù)并代入(6)式， 通過求解線性方程組得到一組PΛne; 第四步， 由離子的布居PΛne計(jì)算等離子體的平均離化度〈Z〉和自由電子密度Ne; 第五步， 根據(jù)新的Ne， 重復(fù)第一步到第四步; 最后， 當(dāng)Ne的變化小于某一數(shù)值(例如千分之一)時(shí)， 就得到了一組收斂的PΛne和

2.3 等離子體效應(yīng)

在稠密等離子體中， 某個離子相鄰的離子和電子產(chǎn)生的電場會對離子的電離能和軌道能產(chǎn)生影響.當(dāng)?shù)入x子體密度足夠高時(shí)， 會產(chǎn)生壓制電離現(xiàn)象.這里， 用Stewart和Pyatt的方法[11]計(jì)算等離子體環(huán)境造成的閾值下降 ΔIne:

這里，IH和a0分別是Rydberg常數(shù)和氫原子半徑;rne=[3(Z-ne)/(4πNe)]1/3是 離 子 球 半 徑;D=[4π(〈Z〉+〈Z2〉)Ni/Te]-1/2是Debye半徑.〈Z〉是等離子體的平均離化度，〈Z2〉是離化度分布平方的平均值.

3 計(jì)算結(jié)果

為了檢驗(yàn)MAICRM模型的有效性， 用該模型計(jì)算了Fe， Xe和Au非平衡穩(wěn)態(tài)等離子體的離化度分布.在本文的計(jì)算中， 每個離化度選用了65個單軌道基， 具體為:nlj(n≤5;l=0，···，n-1) 和nl(5＜n≤10;l=0，···，n-1).考慮閾值下降后，各離化度離子的最大軌道數(shù)nmax≤65.

3.1 Fe等離子體

圖1為計(jì)算的Fe等離子體的平均離化度和離化度分布與其他理論和實(shí)驗(yàn)結(jié)果的比較.

圖1 Fe等離子體的平均離化度和離化度分布的比較.括號里的數(shù)值是實(shí)驗(yàn)的誤差范圍和GALAXY計(jì)算值的變化范圍Fig.1.The comparison of the mean ionization and charge state distribution (CSD) of Fe plasma.The data in the parenthesis are the uncertainties of experimental value and the variation range of GALAXY codes.

圖1 中， 實(shí)驗(yàn)離化度分布[12]通過HULLAC[13]和Kelly[14]計(jì)算的幾千條譜線擬合實(shí)驗(yàn)的吸收譜得到.圖1中， Cloudy[15]， GALAXY[16]和MAICRM計(jì) 算 的 平 均 離 化 度，和在實(shí)驗(yàn)測量誤差范圍內(nèi)，〈Z〉Exp.=16.1±0.2; Flychk[17]和GALAXY計(jì)算的平均離化 度0.2比實(shí)驗(yàn)值略高.這里Cloudy是一種組態(tài)和能級混合的碰撞輻射模型， GALAXY是組態(tài)層次的碰撞輻射模型， Flychk是超組態(tài)層次的碰撞輻射模型.Flychk和MAICRM計(jì)算的離化度分布在電子溫度Te為90—150 eV范圍內(nèi)， 對溫度不敏感.MAICRM計(jì)算的Te=210eV 的離化度分布與實(shí)驗(yàn)符合最好.對于布居最大的Fe15+―Fe18+離子，F(xiàn)lychk和GALAXY計(jì)算的Fe17+和Fe18+離子的布居比MAICRM和Cloudy的高， 但是Fe16+離子的布居要低， 這也是Flychk和Cloudy計(jì)算的平均離化度高的原因.其他三種模型計(jì)算結(jié)果的差異主要來自原子參數(shù)和不同組態(tài)的選擇， 而MAICRM和它們的差別來源于原子參數(shù)和建模方法的差別.

3.2 Xe等離子體

圖2是不同理論方法計(jì)算的Xe等離子體在Te=415eV，Ni=4.75×1018cm-3時(shí)的離化度分布.

圖2 Xe等離子體在 T e=415eV ，Ni=4.75×1018cm-3時(shí)的離化度分布.括號里的數(shù)是實(shí)驗(yàn)的誤差范圍Fig.2.The comparison of the mean ionization and CSD of Xe plasma at T e=415eV ， N i=4.75×1018cm-3.The data in parenthesis is the experimental uncertainties.

圖2 中， 實(shí)驗(yàn)的電子溫度Te=415eV ， 電子密度Ne=(1.30±0.05)×1020cm-3和平均離化度〈Z〉Exp.=27.4±1.5是通 過 擬 合 湯 姆 孫 散 射 譜(Thomson scattering)得到的[18].AVERROES是基于超組態(tài)的理論方法[19].SOSA fit是用自旋-軌道分裂計(jì)算譜線， 并通過擬合實(shí)驗(yàn)譜得到離化度分布[20].三個模型計(jì)算的平均離化度，〈Z〉SOSAfit=26.5 ，〈Z〉A(chǔ)VERROES=26.8 和〈Z〉MAICRM=27.0 ， 都在實(shí)驗(yàn)測量誤差范圍內(nèi)，〈Z〉Exp.=27.4±1.5.相比于另外兩個理論計(jì)算結(jié)果， MAICRM 的離化度分布更偏向高離化度分布.原因是MAICRM用每個離化度所有離子的軌道占據(jù)平均值來計(jì)算電離/復(fù)合速率系數(shù)， 原則上包含了基態(tài)、單激發(fā)態(tài)和多激發(fā)態(tài)的貢獻(xiàn); 而Cloudy和Galaxy 通過挑選每個離化度重要的能態(tài)并計(jì)算其布居， 由于計(jì)算量的限制， 只能挑選有限的激發(fā)態(tài)， 因此考慮的激發(fā)態(tài)，特別是多激發(fā)態(tài)更全面， 導(dǎo)致激發(fā)軌道的平均占據(jù)數(shù)更高， 更容易電離， 因此離化度偏高.

3.3 Au等離子體

圖3是MAICRM方法計(jì)算的Au等離子體在不同狀態(tài)的〈Z〉和離化度分布與實(shí)驗(yàn)測量[21]的比較.

圖3中， 實(shí)驗(yàn)的電子溫度和電子密度是由湯姆孫散射譜和Au+K:Cl混合靶的K殼層的發(fā)射譜診斷得到的.實(shí)驗(yàn)的平均離化度〈Z〉和離化度分布是通過理論計(jì)算的幾百萬條Au41+―Au53+離子的譜線擬合實(shí)驗(yàn)譜得到的.圖3(a)，(b)和(c)中，MAICRM計(jì)算的離化度分布與實(shí)驗(yàn)值符合.由于擬合實(shí)驗(yàn)譜時(shí)沒有考慮離化度小于Au41+的離子，所以圖3(d)中Au37+―Au40+離子的布居沒有對應(yīng)的實(shí)驗(yàn)值.圖3(e)中， 除了{(lán)Te= 2.0 keV，Ne= 6 ×1020cm—3}和{Te= 1.7 keV，Ne= 7 × 1020cm—3}兩個狀態(tài)， MAICRM計(jì)算的〈Z〉都與實(shí)驗(yàn)值在誤差范圍內(nèi)符合.對于上述兩個狀態(tài)， MAICRM計(jì)算的〈Z〉略低于實(shí)驗(yàn)值.

4 討 論

MAICRM模型中有兩個近似， 第一個近似認(rèn)為同一離化度離子的單軌道速率系數(shù)相同; 第二個近似將(4)式拆分為(5)式和(6)式.第一個近似與SCA模型中同一離化度、不同超組態(tài)的單軌道速率系數(shù)相同的近似相類似， 不同的軌道占據(jù)對速率系數(shù)的影響主要反映在占據(jù)數(shù)和空穴數(shù)上.第二個近似將和PΛne的求解分為兩步， 其適用條件是激發(fā)/退激發(fā)過程比電離/復(fù)合過程快， 或者和PΛne的相互影響是一個小量.對于中、高密度等離子體， 由于碰撞激發(fā)/退激發(fā)過程相對電離/復(fù)合過程是一個快過程， 因此， 通過迭代求解(5)式和(6)式可以很快得到收斂的和PΛne.對于低密度且光學(xué)薄的等離子體， 電子碰撞激發(fā)和光激發(fā)速率系數(shù)遠(yuǎn)小于自發(fā)輻射速率系數(shù)， 各離化度離子主要處于基態(tài)占據(jù)， 激發(fā)軌道的占據(jù)數(shù)是一個接近零的小量， 因此電離/復(fù)合過程對軌道占據(jù)的影響不再是一個小量， 這種情況下(4)式需要按不同的方式拆分， 這將另文報(bào)道.

DCA/SCA模型對于不同狀態(tài)的等離子體， 需要選擇不同的組態(tài)/超組態(tài).由于計(jì)算時(shí)間限制，人們通常選擇基態(tài)和最外兩/三個殼層的單激發(fā)、雙激發(fā)態(tài)， 如何保證完備性是該方法的一個難點(diǎn).AA模型用一個平均離子模擬等離子體， 細(xì)致的離化度分布的信息無法獲得.MAICRM模型用多個平均離子模擬等離子體， 可以獲得離化度分布以及不同離化度的軌道占據(jù)數(shù)， 因此相對AA模型更細(xì)致， 而且避免了DCA/SCA選組態(tài)的困難.此外，AA模型只需計(jì)算no個軌道的速率方程; DCA/SCA需要計(jì)算nZ+1個離化度， 每個離化度需要考慮幾百甚至幾千個組態(tài)/超組態(tài)nZ+1，C， 因此計(jì)算量是AA模型的nZ+1，C倍; MAICRM模型考慮nZ+1個離子， 每個離子考慮no個軌道， 因此計(jì)算量約為AA的nZ+1倍， 相比于DCA/SCA模型極大地減小了計(jì)算量.

圖3 MAICRM計(jì)算的不同狀態(tài)下Au等離子體的 〈 Z〉 和離化度分布與實(shí)驗(yàn)測量[21]的比較.(a)， (b)， (c)和(d)中括號里的數(shù)是實(shí)驗(yàn)誤差， (e)中電子溫度Te， 輻射場溫度Tr和電子密度Ne的單位分別是keV， eV和1020 cm—3Fig.3.The comparisons of the mean ionization 〈 Z〉 and CSD of Au plasma between MAICRM and the experimental results[21].In panels (a)， (b)， (c) and (d) the data in parenthesis are the experimental uncertainties.In panel (e) the units of electron temperature Te， radiation temperature Tr and electron density Ne are keV， eV and 1020 cm—3 respectively.

5 結(jié) 論

MAICRM通過迭代求解速率方程組得到平均離子的布居PΛne和軌道占據(jù)數(shù)本文報(bào)道了用該方法計(jì)算的Fe， Xe和Au非平衡等離子體的〈Z〉及離化度分布， 計(jì)算的〈Z〉和離化度分布與實(shí)驗(yàn)測量符合， 精度達(dá)到DCA/SCA的結(jié)果， 且計(jì)算量約為AA模型的nZ+1倍， 遠(yuǎn)低于DCA/SCA的計(jì)算量.因此， MAICRM模型有望成為輻射流體程序所需的在線計(jì)算細(xì)致原子模型.

附錄A 單軌道速率系數(shù)

單軌道速率系數(shù)R是通過截面對電子和光子分布函數(shù)進(jìn)行平均得到的， 正、逆過程速率系數(shù)滿足細(xì)致平衡原理.如果沒有特別說明， 電子分布函數(shù)為Maxwellian分布.注意這里的單軌道速率系數(shù)不包括軌道的占據(jù)數(shù)和空穴數(shù).

A1 自發(fā)輻射速率系數(shù)RSP， 受激輻射速率系數(shù)RST， 光激發(fā)速率系數(shù)RPE

u軌道到l軌道的愛因斯坦自發(fā)輻射速率系數(shù)為

這里，e是電子電荷，c是光速，me是電子質(zhì)量，h是普朗克常數(shù)，分別是軌道間的躍遷能和振子強(qiáng)度.自發(fā)輻射速率系數(shù)、受激輻射速率系數(shù)和光激發(fā)速率系數(shù)分別為

這里，f(ν) 是光子分布函數(shù); 平衡狀態(tài)下光子為普朗克分布， 即f(ν)=1/(ehv/T-1) ;φ(hν) 是譜線的線型因子.

A2 碰撞激發(fā)速率系數(shù)RCE， 碰撞退激發(fā)速率系數(shù)RCD

碰撞激發(fā)截面采用Van Regemorter公式[22]

這里，εf是入射電子能量，εul是激發(fā)能;g(x) 是有效Gaunt因子， 采用Mewe的公式[23]， 即g(x)=A+(B/x)+(C/x2)+Dln(x)，B=C=0 ，D=0.28 ， 對 于Δn=0的躍遷，A=0.15 ， 對于 Δn/=0 的躍遷，A=0.6.碰撞激發(fā)和退激發(fā)速率系數(shù)由碰撞激發(fā)截面與自由電子分布函數(shù)卷積得到:

A3 光電離速率系數(shù)RPI， 輻射復(fù)合速率系數(shù)RRR

光電離截面采用Kramers公式[24]:

這里，ni是電離電子的軌道主量子數(shù);Ine，i是電離閾值， 這里用軌道能近似;hν是光子能量.光電離和輻射復(fù)合速率系數(shù)為

A4 碰撞電離速率系數(shù)RCI， 三體復(fù)合速率系數(shù)R3R

碰撞電離截面用Lotz公式[25]:

這里C值為2.77.碰撞電離和三體復(fù)合速率系數(shù)為

A5 自電離速率系數(shù)RAI， 電子俘獲速率系數(shù)REC

這里Zeff和ni是被俘獲電子的有效核電荷數(shù)和主量子數(shù).俘獲電子的能量寬度用估算.