一个快速模拟热稠密非平衡等离子体的碰撞辐射模型*

韩小英 李凌霄 戴振生 郑无敌 古培俊 吴泽清

(北京应用物理与计算数学研究所， 北京 100094)

在实验室和天体等离子体研究中， 原子/离子的激发、退激发以及电离、复合过程对等离子体的电离和能量平衡有着重要的影响.在激光等离子体作用的辐射流体模拟中， 需要在线计算等离子体的平均离化度和吸收/发射系数.在现有的计算能力下， 通常采用比较简单的平均原子(average atom， AA)模型进行在线计算.随着实验技术和计算能力的发展， 急需发展能够在线计算的细致非平衡原子模型.本文报道了最新发展的多离化度平均离子碰撞辐射模型(multi-average ion collisional-radiative model， MAICRM).该模型用一个平均离子模拟等离子体中某一离化度所有离子的平均轨道占据数和布居， 即每个平均离子的轨道占据数为该离化度所有离子的轨道占据数的平均; 平均离子的布居等于该离化度离子的布居和.平均离子的轨道占据数和布居通过迭代求解速率方法得到.用该模型计算了Fe， Xe和Au非平衡等离子体的离化度分布， 计算结果与细致组态和超组态模型以及实验测量符合， 而计算量相对于细致组态/超组态大大降低.预期该方法能与辐射流体程序耦合， 实现细致非平衡原子模型的在线计算.

1 引 言

在实验室和天体等离子体的辐射流体模拟中，等离子体的离化度分布是其中关键的物理量[1，2].等离子体的电离平衡和离化度分布主要由原子(离子)的激发/退激发和电离/复合过程决定.当密度高时， 等离子体接近局域热动平衡状态(local thermodynamic equilibrium， LTE)， 离子的轨道占据数和布居满足Saha-Boltzmann分布; 当密度低时， 等离子体处于冕区平衡状态(coronal equilibrium， CE); 对于中等密度的等离子体， 则需要通过求解速率方程得到各离子的非平衡分布， 通常称为碰撞-辐射(collisional-radiative， CR)模型.

目前模拟非平衡等离子体的方法主要有以下几种(从考虑的能态由粗到细): 平均原子模型(average atom， AA)[3]、超组态方法(super-configuration accounting， SCA)[4]、细致组态方法(detailed configuration accounting， DCA)[5]、细致能级方法(detailed level accounting， DLA)[6].其中AA模型用一个平均原子模拟等离子体中各离子轨道的平均占据数给出平均离化度.AA模型的计算时间最短， 但是无法给出细致的离化度分布.DCA，SCA和DLA方法能给出大量能态(几千或者几万个能态)的细致分布， 但是计算时间比AA模型高几个量级.此外， 上述细致原子模型对于不同温度密度下的等离子体， 通常只能考虑有限的能态， 如何保证所选能态的完备性一直是该领域的一个难题.

在激光等离子体作用中， 等离子体要经历从低温高密到高温低密的宽范围状态， 目前通常采用AA模型在线计算等离子体的平均离化度和吸收、发射系数.从吸收和发射总量来看， AA模型与实验测量符合[7，8].近年来， 随着实验技术的发展， 人们不仅能测量等离子体的发射总量， 而且能够进行谱分辨的细致测量; 理论上也发展了所谓的“DCA”在线非平衡原子模型方法， 但是这里的“DCA”是用主量子数n刻画组态[9]， 类似于n层次的超组态方法.文献[10]比较了不同版“DCA”模型计算的发射总量， 发现计算结果有倍数差别.

针对细致、快速在线非平衡原子模型的需求，本研究组参考AA模型的物理图像， 基于细致组态的速率方程， 发展了一种新的非平衡原子碰撞辐射模型， 即多离化度平均离子碰撞辐射模型(multiaverage ion collisional-radiative model， MAICRM).首先， 用一个平均离子来模拟某一离化度所有组态的轨道平均占据数， 该占据数通过求解轨道之间的激发和退激发速率方程得到.然后， 基于平均离子的轨道占据数， 通过求解离子之间的电离和复合速率方程， 得到各平均离子的布居.有了平均离子的轨道占据数和布居， 就可以得到等离子体的平均离化度和离化度分布以及细致的吸收和发射谱.该方法采用了两个近似: 第一个近似认为同一离化度离子的单轨道速率系数(不包含轨道的占据数和空穴数)相等; 第二个近似将束缚轨道之间的跃迁(激发和退激发过程)和束缚轨道与连续轨道之间的跃迁(电离和复合过程)分两步计算.

在DCA/SCA模型中， 每个离化度需要考虑几百甚至几千个组态/超组态; AA模型只需考虑一个离化度的十个或者几十个轨道之间的跃迁;MAICRM模型中， 最多只需考虑Z+1 个离子之间的电离、复合过程， 每个离子内计算nmax个轨道之间的激发、退激发过程， 因此MAICRM相对于DCA/SCA模型能极大地减少计算量， 其计算量约为AA模型的Z+1 倍.

本文报道了用该方法计算的Fe， Xe和Au非平衡等离子体的离化度分布和平均离化度， 计算结果与DCA/SCA计算结果以及实验测量符合， 因此该方法有望在未来与辐射流体耦合， 实现在线细致非平衡原子模型计算.

2 理论方法

众所周知， AA模型是用一个平均原子模拟所有离子轨道的平均占据， 该模型计算的平均离化度(核电荷数-束缚轨道占据数)、总发射量等与实验测量基本相符， 但是无法获得更细致的离化度分布等信息.从AA模型图像出发， 考虑用一个平均离子模拟某个离化度所有离子的轨道平均占据数，从而可以得到更细致的离化度分布和不同离化度轨道的平均占据数.基于该物理图像， 可以从细致组态的速率方程出发， 通过对同一离化度所有组态的轨道占据数求平均， 得到各离化度离子的轨道平均占据数和布居的速率方程.

2.1 细致组态的速率方程

有ne个束缚电子的组态记为表示该组态各轨道的电子占据数， 组态的布居为PKne;gi和εi分别表示i轨道的统计权重和轨道能.该组态布居的速率方程如下:

式中， RE/D为激发/退激发速率系数， 包括电子碰撞和光过程.RI包括电子碰撞电离和光电离速率系数， RR包括辐射复合和三体复合速率系数，RAI为自电离速率系数， REC为电子俘获速率系数.各速率系数的计算公式详见附录A.注意， 这里的速率系数 R 为单轨道速率系数， 不包含涉及轨道的占据数和空穴数.

2.2 平均离子的速率方程

将(1)式两边乘以Kne组态的n轨道的占据数并对该离化度所有组态求和， 得到

求 和 过 程 中， 对 于u→l的 跃 迁， 利 用 了 关 系 式这里认为同一离化度各组态的单电子轨道速率系数 R 相同， 即 R 只与离化度有关.该近似与SCA方法中认为同一超组态的单轨道速率系数相同类似.

进一步， 定义

其中，PΛne是所有具有ne个束缚电子的组态(Kne)的布居之和;是同一离化度所有组态(Kne)n轨道上占据数的平均值.据此， 定义一个平均离子它由以下三个物理量刻画: 1)有ne个束缚电子; 2)布 居 为PΛne; 3)n轨 道 上 占 据 数 为为Λne离 子 的n轨 道 上 的 空 穴 数.为Λne离子u和l轨道上占据数之积.将(3)式代入(2)式可得:

(4)式即为平均离子Λne的布居PΛne和轨道占据数的速率方程.在(4)式中是耦合的，很难直接通过(4)式求解考虑到一般况下， 电子的束缚-束缚跃迁速率比束缚-连续跃迁速率快， 轨道上电子占据数通常由激发和退激发过程决定， 而离化度的布居通常由电离和复合过程决定.因此， 将(4)式拆为下面两个式子:(5)式为离子轨道占据数的速率方程， 包括轨道之间的激发和退激发过程.(6)式为离子布居的速率方程， 包括电离和复合过程.通过迭代求解(5)式和(6)式， 可以得到一组收敛的具体的计算过程为: 第一步， 根据给定的等离子体参数， 包括电子温度(Te)、辐射场温度(Tr)、离子温度(Ti)和电子密度(Ne)， 计算激发/退激发速率系数 RE/D; 第二步， 将激发/退激发速率系数代入(5)式， 假设每个离子初始占据数为基态占据， 再用隐式牛顿迭代法得到每个离子收敛的第三步， 根据每个离子的轨道占据数计算电离/复合速率系数并代入(6)式， 通过求解线性方程组得到一组PΛne; 第四步， 由离子的布居PΛne计算等离子体的平均离化度〈Z〉和自由电子密度Ne; 第五步， 根据新的Ne， 重复第一步到第四步; 最后， 当Ne的变化小于某一数值(例如千分之一)时， 就得到了一组收敛的PΛne和

2.3 等离子体效应

在稠密等离子体中， 某个离子相邻的离子和电子产生的电场会对离子的电离能和轨道能产生影响.当等离子体密度足够高时， 会产生压制电离现象.这里， 用Stewart和Pyatt的方法[11]计算等离子体环境造成的阈值下降 ΔIne:

这里，IH和a0分别是Rydberg常数和氢原子半径;rne=[3(Z-ne)/(4πNe)]1/3是 离 子 球 半 径;D=[4π(〈Z〉+〈Z2〉)Ni/Te]-1/2是Debye半径.〈Z〉是等离子体的平均离化度，〈Z2〉是离化度分布平方的平均值.

3 计算结果

为了检验MAICRM模型的有效性， 用该模型计算了Fe， Xe和Au非平衡稳态等离子体的离化度分布.在本文的计算中， 每个离化度选用了65个单轨道基， 具体为:nlj(n≤5;l=0，···，n-1) 和nl(5＜n≤10;l=0，···，n-1).考虑阈值下降后，各离化度离子的最大轨道数nmax≤65.

3.1 Fe等离子体

图1为计算的Fe等离子体的平均离化度和离化度分布与其他理论和实验结果的比较.

图1 Fe等离子体的平均离化度和离化度分布的比较.括号里的数值是实验的误差范围和GALAXY计算值的变化范围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 中， 实验离化度分布[12]通过HULLAC[13]和Kelly[14]计算的几千条谱线拟合实验的吸收谱得到.图1中， Cloudy[15]， GALAXY[16]和MAICRM计 算 的 平 均 离 化 度，和在实验测量误差范围内，〈Z〉Exp.=16.1±0.2; Flychk[17]和GALAXY计算的平均离化 度0.2比实验值略高.这里Cloudy是一种组态和能级混合的碰撞辐射模型， GALAXY是组态层次的碰撞辐射模型， Flychk是超组态层次的碰撞辐射模型.Flychk和MAICRM计算的离化度分布在电子温度Te为90—150 eV范围内， 对温度不敏感.MAICRM计算的Te=210eV 的离化度分布与实验符合最好.对于布居最大的Fe15+―Fe18+离子，Flychk和GALAXY计算的Fe17+和Fe18+离子的布居比MAICRM和Cloudy的高， 但是Fe16+离子的布居要低， 这也是Flychk和Cloudy计算的平均离化度高的原因.其他三种模型计算结果的差异主要来自原子参数和不同组态的选择， 而MAICRM和它们的差别来源于原子参数和建模方法的差别.

3.2 Xe等离子体

图2是不同理论方法计算的Xe等离子体在Te=415eV，Ni=4.75×1018cm-3时的离化度分布.

图2 Xe等离子体在 T e=415eV ，Ni=4.75×1018cm-3时的离化度分布.括号里的数是实验的误差范围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 中， 实验的电子温度Te=415eV ， 电子密度Ne=(1.30±0.05)×1020cm-3和平均离化度〈Z〉Exp.=27.4±1.5是通 过 拟 合 汤 姆 孙 散 射 谱(Thomson scattering)得到的[18].AVERROES是基于超组态的理论方法[19].SOSA fit是用自旋-轨道分裂计算谱线， 并通过拟合实验谱得到离化度分布[20].三个模型计算的平均离化度，〈Z〉SOSAfit=26.5 ，〈Z〉AVERROES=26.8 和〈Z〉MAICRM=27.0 ， 都在实验测量误差范围内，〈Z〉Exp.=27.4±1.5.相比于另外两个理论计算结果， MAICRM 的离化度分布更偏向高离化度分布.原因是MAICRM用每个离化度所有离子的轨道占据平均值来计算电离/复合速率系数， 原则上包含了基态、单激发态和多激发态的贡献; 而Cloudy和Galaxy 通过挑选每个离化度重要的能态并计算其布居， 由于计算量的限制， 只能挑选有限的激发态， 因此考虑的激发态，特别是多激发态更全面， 导致激发轨道的平均占据数更高， 更容易电离， 因此离化度偏高.

3.3 Au等离子体

图3是MAICRM方法计算的Au等离子体在不同状态的〈Z〉和离化度分布与实验测量[21]的比较.

图3中， 实验的电子温度和电子密度是由汤姆孙散射谱和Au+K:Cl混合靶的K壳层的发射谱诊断得到的.实验的平均离化度〈Z〉和离化度分布是通过理论计算的几百万条Au41+―Au53+离子的谱线拟合实验谱得到的.图3(a)，(b)和(c)中，MAICRM计算的离化度分布与实验值符合.由于拟合实验谱时没有考虑离化度小于Au41+的离子，所以图3(d)中Au37+―Au40+离子的布居没有对应的实验值.图3(e)中， 除了{Te= 2.0 keV，Ne= 6 ×1020cm—3}和{Te= 1.7 keV，Ne= 7 × 1020cm—3}两个状态， MAICRM计算的〈Z〉都与实验值在误差范围内符合.对于上述两个状态， MAICRM计算的〈Z〉略低于实验值.

4 讨 论

MAICRM模型中有两个近似， 第一个近似认为同一离化度离子的单轨道速率系数相同; 第二个近似将(4)式拆分为(5)式和(6)式.第一个近似与SCA模型中同一离化度、不同超组态的单轨道速率系数相同的近似相类似， 不同的轨道占据对速率系数的影响主要反映在占据数和空穴数上.第二个近似将和PΛne的求解分为两步， 其适用条件是激发/退激发过程比电离/复合过程快， 或者和PΛne的相互影响是一个小量.对于中、高密度等离子体， 由于碰撞激发/退激发过程相对电离/复合过程是一个快过程， 因此， 通过迭代求解(5)式和(6)式可以很快得到收敛的和PΛne.对于低密度且光学薄的等离子体， 电子碰撞激发和光激发速率系数远小于自发辐射速率系数， 各离化度离子主要处于基态占据， 激发轨道的占据数是一个接近零的小量， 因此电离/复合过程对轨道占据的影响不再是一个小量， 这种情况下(4)式需要按不同的方式拆分， 这将另文报道.

DCA/SCA模型对于不同状态的等离子体， 需要选择不同的组态/超组态.由于计算时间限制，人们通常选择基态和最外两/三个壳层的单激发、双激发态， 如何保证完备性是该方法的一个难点.AA模型用一个平均离子模拟等离子体， 细致的离化度分布的信息无法获得.MAICRM模型用多个平均离子模拟等离子体， 可以获得离化度分布以及不同离化度的轨道占据数， 因此相对AA模型更细致， 而且避免了DCA/SCA选组态的困难.此外，AA模型只需计算no个轨道的速率方程; DCA/SCA需要计算nZ+1个离化度， 每个离化度需要考虑几百甚至几千个组态/超组态nZ+1，C， 因此计算量是AA模型的nZ+1，C倍; MAICRM模型考虑nZ+1个离子， 每个离子考虑no个轨道， 因此计算量约为AA的nZ+1倍， 相比于DCA/SCA模型极大地减小了计算量.

图3 MAICRM计算的不同状态下Au等离子体的 〈 Z〉 和离化度分布与实验测量[21]的比较.(a)， (b)， (c)和(d)中括号里的数是实验误差， (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 结 论

MAICRM通过迭代求解速率方程组得到平均离子的布居PΛne和轨道占据数本文报道了用该方法计算的Fe， Xe和Au非平衡等离子体的〈Z〉及离化度分布， 计算的〈Z〉和离化度分布与实验测量符合， 精度达到DCA/SCA的结果， 且计算量约为AA模型的nZ+1倍， 远低于DCA/SCA的计算量.因此， MAICRM模型有望成为辐射流体程序所需的在线计算细致原子模型.

附录A 单轨道速率系数

单轨道速率系数R是通过截面对电子和光子分布函数进行平均得到的， 正、逆过程速率系数满足细致平衡原理.如果没有特别说明， 电子分布函数为Maxwellian分布.注意这里的单轨道速率系数不包括轨道的占据数和空穴数.

A1 自发辐射速率系数RSP， 受激辐射速率系数RST， 光激发速率系数RPE

u轨道到l轨道的爱因斯坦自发辐射速率系数为

这里，e是电子电荷，c是光速，me是电子质量，h是普朗克常数，分别是轨道间的跃迁能和振子强度.自发辐射速率系数、受激辐射速率系数和光激发速率系数分别为

这里，f(ν) 是光子分布函数; 平衡状态下光子为普朗克分布， 即f(ν)=1/(ehv/T-1) ;φ(hν) 是谱线的线型因子.

A2 碰撞激发速率系数RCE， 碰撞退激发速率系数RCD

碰撞激发截面采用Van Regemorter公式[22]

这里，εf是入射电子能量，εul是激发能;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.碰撞激发和退激发速率系数由碰撞激发截面与自由电子分布函数卷积得到:

A3 光电离速率系数RPI， 辐射复合速率系数RRR

光电离截面采用Kramers公式[24]:

这里，ni是电离电子的轨道主量子数;Ine，i是电离阈值， 这里用轨道能近似;hν是光子能量.光电离和辐射复合速率系数为

A4 碰撞电离速率系数RCI， 三体复合速率系数R3R

碰撞电离截面用Lotz公式[25]:

这里C值为2.77.碰撞电离和三体复合速率系数为

A5 自电离速率系数RAI， 电子俘获速率系数REC

这里Zeff和ni是被俘获电子的有效核电荷数和主量子数.俘获电子的能量宽度用估算.