李树

(北京应用物理与计算数学研究所,北京 100094)(2018年5月10日收到;2018年8月16日收到修改稿)

高温全电离等离子体的辐射输运问题中,光子与电子的Compton散射与逆Compton散射是其中重要的特性,光子与相对论麦克斯韦电子散射的描述及截面的计算非常复杂且费时.本文提出了一种用于模拟计算光子与相对论麦克斯韦速度分布电子散射截面的蒙特卡罗计算方法.给出了各步骤的具体实现办法,推导了对应的计算公式,研究了相对论电子速率抽样方法,编写了光子与相对论电子散射的微观截面的蒙特卡罗计算程序.开展了高温全电离等离子体中,不同能量光子与不同温度电子散射的微观散射截面计算和分析.模拟计算结果显示,在电子温度低于25 keV情况下,本文方法与多重数值积分方法的计算结果非常接近;但随着电子温度继续升高,二者差异逐渐增大并较明显,经分析,可能是本文方法目前的电子速率抽样偏差所致,希望将来能够找到更好的相对论电子速率抽样方法以克服此缺陷.

1 引 言

恒星内部、惯性约束核聚变(inertial confinement fusion,ICF)聚变区、热核装置等高温等离子体系统中的温度高达上亿摄氏度.在如此高温度的环境下,轻物质完全电离,由热辐射产生的光子能量范围主要在103eV量级.对于完全电离的等离子体,能量为103eV量级的光子与物质(电子)之间的能量交换过程包括:轫致辐射、逆轫致吸收、Compton散射和逆Compton散射.于敏院士[1]在“高温等离子体中物质与光之间的能量传递过程”一文中探讨并得出结论“电子能量损失主要通过(逆)Compton散射”.因此,在研究恒星或ICF中的辐射输运问题时,Compton及逆Compton散射问题是需要重点关注的对象.

当上述等离子体系统处于完全热平衡态(电子、离子和辐射的温度相同)时,电子温度可达10 keV以上;某些情况下,系统处于非平衡态(电子、离子、辐射具有各自的温度),电子温度可能达到近100 keV.在此非平衡态情况下,光子与电子的相互作用不仅需考虑电子的热运动速度,而且必需在相对论体系下考虑光子与电子的相互作用[2].

利用蒙特卡罗方法模拟高温等离子体中的辐射输运问题时,需要知道任意能量的光子与任意温度电子的相互作用(吸收、散射)截面,以及发生散射后出射光子的能量和角度.本文讨论光子与电子的散射截面计算问题,散射后光子和电子的能谱、角度谱计算将在后续的研究中讨论.

光子与静止自由电子发生Compton散射的微分截面Klein-Nishina公式[3]为式中re为电子经典半径;ε为以电子静止能量为单位的入射光子能量,即ε=hν/(m0c2),其中hν为入射光子能量,m0为静止电子质量,c为光速(约300 cm/sh,sh是本文的时间单位,1 sh=10−8s);ϑ为入射光子与出射光子夹角(偏转角);Ω为立体角.

对(1)式关于立体角积分,得到能量为hν的光子与1个静止自由电子发生散射的概率(微观散射总截面)公式[4]

当ε较小时,(2)式的计算容易产生数值误差,且该式中的对数计算较费时.故通常采用如下的Hastings[5]近似公式代替(2)式计算微观散射总截面:

式中η= 1+0.222037ε.这 里 称(3)式 为KNHastings公式.

在特别设定的坐标系(入射光子沿Z轴飞行、出射光子位于XZ平面)下,光子与速度分布函数为f(u)的各向同性电子发生散射的微分截面计算公式[2]为式中ν和ν′分别为入射、出射光子频率;ζ=cosϑ=Ω ·Ω′为光子偏转角余弦;u为电子速率;为电子的Lorentz因子;µ=cosα为与电子运动方向相关的第1个积分角度变量,其中α为电子运动方向与入射光子方向(Z轴)的夹

角;ε=hν/(m0c2);ε′=hν′/(m0c2);函数g(u,µ)定义为

ϕ为电子运动的方位角.δ函数是为了满足散射的动量守恒引入的,且仅当

时,(5)式中的δ可以消除掉.电子速度分布f(u)可取任意概率密度函数,这里取相对论麦克斯韦-玻尔兹曼(relativistic Maxwellian-Boltzman,RMB)速度分布函数(详见第3节).

关于(4)式的计算比较复杂,已有多篇文献介绍如何计算该式.这些文献[6−15]的基本思路大多是:首先将(4)式中的积分核函数展开(幂级数、勒让得函数),然后采用多重数值积分的方法计算散射截面以及在此基础上计算确定论输运模拟方法[16]所需的分群能谱与角度谱、或者蒙特卡罗输运模拟所需的散射后出射光子的能量与方向[17,18].

本文不讨论数值计算(4)式,而是利用蒙特卡罗方法,首先逐个直接模拟光子与相对论电子的相互作用,得到若干光子与单个电子散射的概率值,然后取这些概率值的统计平均值,最后得到一定频率的光子与一定温度的电子发生散射的概率(截面).因此,本文方法的实质是利用蒙特卡罗方法计算多维复杂核函数的积分问题(注:计算散射截面为6维,计算散射后的能谱角度谱为8维).

2 光子与相对论电子散射的蒙特卡罗方法

图1 实验室系光子与温度为Te的电子散射示意图Fig.1.Photon scattering with electrons at temperature Tein lab coordinate.

问题描述:如图1所示,设实验室坐标系X0Y0Z0(简称0系)中,t时刻有一光子(hν0,Ω0)(光子能量和飞行方向),求其与r0处温度为Te的电子发生散射的概率(微观散射截面)及散射后光子和电子的能量、角度分布.由于篇幅的原因,本文只讨论散射概率的计算方法,散射后光子和电子的能量、角度分布计算方法将另文讨论.

蒙特卡罗求解步骤:

步骤1根据电子温度Te抽样单个电子状态(u0e,Ω0e)(电子速率和飞行方向).其中,电子飞行方向Ω0e(Ωx0e,Ωy0e,Ωz0e)为各向同性抽样:

式中ξ1和ξ2为[0,1]之间的随机数.

由此可以求得光子与电子的飞行方向夹角余弦cosα0=Ω0·Ω0e.

电子速率u0e从RMB速率分布函数中抽样得到,后面将针对此问题单独展开讨论.

入射光子与抽样出的单个电子在实验室系中的散射前状态如图2所示.

图2 实验室系中光子与抽样出的单个电子散射前示意图Fig.2.Photon scattering with the sampled electron in lab coordinate.

步骤2设置新笛卡儿坐标系X1Y1Z1(简称1系),其坐标原点与r0重合,X1轴正方向与电子运动方向Ω0e一致,Z1轴垂直于电子运动方向与光子飞行方向所组成的平面.由于0系与1系只存在平移和旋转变换,因此在1系中,光子与X1轴(电子飞行方向)夹角α1=α0,光子在1系中的(碰撞前)状态为

电子在1系中的(碰撞前)状态为

入射光子与单个电子在坐标系1中的散射前状态如图3所示.

图3 坐标系1中光子电子散射前示意图Fig.3.Photon scattering with electron in coordinate 1.

步骤3建立一个随电子运动的笛卡儿坐标系X2Y2Z2(简称2系).坐标系2与坐标系1在t时刻重合,运动方向为X1轴正方向,运动速度为u1e(图4).故电子在2系中始终静止(u2e=0)且位于坐标原点.

图4 坐标系1与坐标系2关系示意图Fig.4.Relationship between coordinate 1 and coordinate 2.

现求光子在2系中的状态.根据狭义相对论,对于光子,(Px,Py,Pz,iE/c)构成一个四度向量,又因光子动量光子能量E=hν,所以,(νΩx,νΩy,νΩz,iν)亦构成四度向量.1系中的四度向量(ν1Ωx1,ν1Ωy1,ν1Ωz1,iν1)与运动坐标系2系中的四度向量(ν2Ωx2,ν2Ωy2,ν2Ωz2,iν2)可以通过Lorentz变换建立关系[19]:式中λ=[1−(u1e/c)2]−1/2.

求解(10)式后得到

代入步骤2已求得的u1e,Ωx1,Ωy1,Ωz1值后得到

由此计算得到光子在2系的频率ν2及飞行方向Ω2(Ωx2,Ωy2,Ωz2),(12)式中α2是光子在2系中的飞行方向与X2轴的夹角(如图5所示,Z1,Z2轴均垂直于纸面).

图5 光子在1系和2系中的飞行方向与坐标轴夹角示意图Fig.5.Photon’s flight directions and angles in coordinate 1 and coordinate 2.

步骤4建立笛卡儿坐标系X3Y3Z3(简称3系).Z3轴与Z2轴相同,X3轴由X2轴以Z2轴为轴旋转角度α2得到,即X3轴与光子在2系中的飞行方向Ω2重合,2系与3系的关系及光子飞行方向如图6所示,其中的Z轴重合且均垂直于纸面.由于2系与3系之间仅存在旋转变换,故光子频率不变,因此光子在3系中的状态为:hν3=hν2,Ωx3=1,Ωy3=0,Ωz3=0.

至此,在3系中,光子与电子的相互作用图像是:一个能量为hν3的光子沿X3轴正方向飞行,在坐标原点处与一个自由静止的电子发生散射,如图7所示.故可以方便地采用Klein-Nishina公式处理,这里的微观散射总截面计算采用(3)式,式中的ε=hν3/(m0c2), 计算出的σs(ε)就是光子(hν0,Ω0)与相对论电子(u0e,Ω0e)发生散射的概率.

图6 坐标系2与坐标系3关系示意图Fig.6.Relationship between coordinate 2 and coordinate 3.

图7 坐标系3中光子与电子散射示意图Fig.7.Photon scattering with electron in coordinate 3.

上述4个步骤描述了入射光子与单个特定(抽样出的)相对论电子的微观散射截面计算方法.由于相对论效应引起的空间收缩,与光子飞行方向相反的电子密度比与光子飞行方向相同的电子密度高,因此,需要将截面σs(ε)乘以空间尺度变化引入的Lorentz因子[20].对于光子与温度为Te的电子发生散射的概率,重复步骤1—步骤4若干(N)次,便可以得到N次平均的结果:

式中σs(ε(n))为第n抽样计算得到的微观散射截面;γ(n)为第n次抽样计算的空间变化Lorentz因子:

式中cosα0=Ω0·Ω0e为光子与电子的飞行方向夹角余弦.

抽样足够多的样本计算后,⟨σs(hν0,Te)⟩将逐渐收敛到某一固定值,该值即可认为是光子(hν0,Ω0)与温度为Te的电子发生散射的概率(微观散射截面).

分析上面的计算流程可以看出,对于散射截面,本文方法的计算精度取决于两个因素:电子速率分布抽样的准确性;抽样的样本数N(收敛性).

3 关于电子速率的抽样

经典麦克斯韦-玻尔兹曼(classical Maxwellian-Boltzman,CMB)速率分布的概率密度函数为

此函数的一种比较新且效率和精度均较高的抽样方法可以参考文献[21].

高温全电离等离子体中,当电子达到热力学平衡态时,温度为Te的电子速率u服从RMB速率分布,其概率密度函数[19,22]为

式中k为玻尔兹曼常数为电子的Lorentz因子;K2(x)为二阶变型贝塞尔函数:

求(16)式的积分非常困难,因此抽样计算速率u的值很难实现.笔者也未能查阅到关于(16)式抽样方法的相关文献.因此,为了解决此问题,本文推导了一种近似RMB(similar relativistic Maxwellian-Boltzman,SRMB)速率分布函数:

式中β=λ−1.关于分布函数的推导及抽样方法不是本文的重点,这里暂且不讨论.

本文利用(18)式作为电子速率分布的抽样函数分析抽样结果.

首先,分析需要多少样本才能使得抽样的电子速度分布收敛.针对(18)式编制抽样程序,计算并统计不同抽样样本条件下的电子速率分布.图8给出了在电子温度为25 keV情况下,样本数分别为104,105和106条件下,统计平均后得到的电子速率概率密度分布及与(18)式解析计算结果的比较情况.可以看出当抽样样本数达到106时,抽样统计结果与解析结果基本一致,由此可以认为结果已收敛.为了尽量排除抽样样本引入的误差,后面的计算结果都是样本数为107的情况下得到的.

图8 三种样本数条件下的电子速率分布及与解析结果比较(Te=25 keV)Fig.8.Electron speed distributions and analysis results with three kinds of sampled particles(Te=25 keV).

然后,分析近似(18)式抽样的电子速率分布的准确性.这里的准确性,指的是与RMB速率分布函数(16)式的符合程度,若符合得越好,则准确性越高.图9给出5种电子温度(1,10,25,50和100 keV)情况下,利用近似(18)式抽样计算得到电子速率分布(SRMB sampled),将其与(16)式的解析结果(RMB analysis)进行比较,同时图中还给出了根据(15)式抽样计算的电子速率分布(CMB sampled).

1)电子温度低于10 keV情况(图9(a)和图9(b))下,相对论效应不明显,利用CMB分布函数、SRMB分布函数抽样计算得到的结果与RMB分布函数的解析结果基本一致.

2)电子温度升高至25 keV情况(图9(c))下,由于相对论效应,利用CMB分布函数抽样计算得到的结果与RMB分布函数的解析结果有较明显的差异;然而,利用SRMB分布函数抽样计算得到的结果与RMB分布函数的解析结果仍有较好的符合度.

图9 不同电子温度下电子速率分布 (a)1 keV;(b)10 keV;(c)25 keV;(d)50 keV;(e)100 keVFig.9.Electron speed distributions under different electron temperatures:(a)1 keV;(b)10 keV;(c)25 keV;(d)50 keV;(e)100 keV.

3)随着电子温度的继续上升 (图9(d)和图9(e)),CMB分布函数抽样计算得到的结果与RMB分布函数的解析结果差异继续增大且非常显著;利用SRMB分布函数抽样计算得到的结果与RMB分布函数的解析结果有差异但小于CMB抽样结果.另外,从图9(d)和图9(e)可以看出,CMB分布抽样出的电子速度有显著超光速(c=300 cm/sh)不合理现象出现,然而,SRMB分布抽样出的电子速度不会出现超光速现象.

综合以上计算结果及分析,可得:样本数106以上,抽样计算结果基本收敛;电子温度不超过25 keV,本文的电子速率分布抽样较准确,电子温度高于25 keV,偏差增大,这可能会影响到散射截面的计算精度.

4 散射截面的数值模拟结果与分析

根据本文第2节介绍的方法,研制了光子-相对论麦克斯韦电子散射模拟程序,该程序可以计算任意光子(能量范围为0—1 MeV)与任意温度电子(温度范围为0—100 keV)的微观散射截面及散射角度谱、能谱.表1是文献[11]中给出的利用数值积分方法计算的光子-电子散射截面.为了分析电子热运动对散射截面影响,表2列出了利用(3)式计算的光子与静止自由电子的散射截面;表3是利用本文方法计算的光子-电子散射截面.为了便于比较,表2和表3在给出散射截面的同时还列出了其与表1中对应值的相对误差其中σtable1,σtable2和σtable3分别是表1、表2和

首先,分析表1可以看出:光子能量越高,截面越小;电子热运动使得截面变小,且电子温度越高、热运动平均速度越大,截面减小越多.

表1 光子-相对论麦克斯韦电子散射截面(文献[11]数值积分计算)Table 1.Photon-Maxwellian electron cross-sections(numerical integration method,Ref.[11]).

表2 光子-静止自由电子散射截面(KN-Hastings近似(3)式计算)Table 2.Photon-still free electron at rest cross-sections(KN-Hastings method,using Eq.(3)calculated).

表3 光子-相对论电子散射截面(本文蒙特卡罗方法计算)Table 3.Photon-Maxwellian electron cross-sections(proposed Monte Carlo method).

表2的截面计算因未考虑电子热运动,故相对于表1其值偏大,且电子温度越高,截面偏大越多;光子能量越高,截面偏大越多;例如:在光子能量为1 MeV、电子温度为100 keV情况下,相对误差达18%左右.

表3本文方法计算的散射截面,在电子温度不高于25 keV情况下与表1的结果很接近,相对误差不超过2%;但是随着电子温度继续升高,偏差逐渐增大,例如:在光子能量为1 MeV、电子温度为100 keV情况下,相对误差也达到7.8%.

分析本文方法在电子温度较高(>25 keV)情况下与国外的数值积分方法差异的原因:前面已经提到过,电子平均速度越大截面越小.然而,从前面的图9(d)和图9(e)可知,在电子温度较高时,本文目前关于电子速率的抽样方法抽出的电子速率整体偏小,这可能就是导致本文方法计算的散射截面偏大的原因.因此,本研究下一步将围绕相对论电子速率的准确抽样目标开展.

5 结 论

本文提出了一种计算光子与相对论麦克斯韦分布电子散射截面的方法.该方法的实质是利用蒙特卡罗方法计算多维复杂核函数的积分.本文方法具有如下特点:第一,在给定的光子能量范围(0—1 MeV)和电子温度范围(0—100 keV)内,不需要针对光子能量和电子温度分群离散,即可以根据单一光子能量和单一电子温度直接计算出对应的散射截面;第二,散射后光子与电子的能量和角度是连续的,可以方便地利用蒙特卡罗方法确定(另文介绍),因此,本文方法特别适用于利用蒙特卡罗方法模拟光子在高温等离子体中的Compton散射和逆Compton散射问题;第三,与数值积分方法计算散射截面相比,本文方法无需做复杂繁琐的数学处理,因此更加直观、简洁.

本文方法在电子温度不超过25 keV情况下的计算精度较高(注:此情况能够满足大多数ICF、热核装置等问题所需);在电子温度更高情况下,可能是目前的电子速率抽样有偏差的原因,导致目前的计算结果与国外多重数值积分结果存在较明显差异;待将来找到更好的相对论电子速率抽样方法后,此缺陷应该能够被克服.