李樹

(北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所,北京 100094)(2018年5月10日收到;2018年8月16日收到修改稿)

高溫全電離等離子體的輻射輸運(yùn)問題中,光子與電子的Compton散射與逆Compton散射是其中重要的特性,光子與相對論麥克斯韋電子散射的描述及截面的計(jì)算非常復(fù)雜且費(fèi)時.本文提出了一種用于模擬計(jì)算光子與相對論麥克斯韋速度分布電子散射截面的蒙特卡羅計(jì)算方法.給出了各步驟的具體實(shí)現(xiàn)辦法,推導(dǎo)了對應(yīng)的計(jì)算公式,研究了相對論電子速率抽樣方法,編寫了光子與相對論電子散射的微觀截面的蒙特卡羅計(jì)算程序.開展了高溫全電離等離子體中,不同能量光子與不同溫度電子散射的微觀散射截面計(jì)算和分析.模擬計(jì)算結(jié)果顯示,在電子溫度低于25 keV情況下,本文方法與多重?cái)?shù)值積分方法的計(jì)算結(jié)果非常接近;但隨著電子溫度繼續(xù)升高,二者差異逐漸增大并較明顯,經(jīng)分析,可能是本文方法目前的電子速率抽樣偏差所致,希望將來能夠找到更好的相對論電子速率抽樣方法以克服此缺陷.

1 引 言

恒星內(nèi)部、慣性約束核聚變(inertial confinement fusion,ICF)聚變區(qū)、熱核裝置等高溫等離子體系統(tǒng)中的溫度高達(dá)上億攝氏度.在如此高溫度的環(huán)境下,輕物質(zhì)完全電離,由熱輻射產(chǎn)生的光子能量范圍主要在103eV量級.對于完全電離的等離子體,能量為103eV量級的光子與物質(zhì)(電子)之間的能量交換過程包括:軔致輻射、逆軔致吸收、Compton散射和逆Compton散射.于敏院士[1]在“高溫等離子體中物質(zhì)與光之間的能量傳遞過程”一文中探討并得出結(jié)論“電子能量損失主要通過(逆)Compton散射”.因此,在研究恒星或ICF中的輻射輸運(yùn)問題時,Compton及逆Compton散射問題是需要重點(diǎn)關(guān)注的對象.

當(dāng)上述等離子體系統(tǒng)處于完全熱平衡態(tài)(電子、離子和輻射的溫度相同)時,電子溫度可達(dá)10 keV以上;某些情況下,系統(tǒng)處于非平衡態(tài)(電子、離子、輻射具有各自的溫度),電子溫度可能達(dá)到近100 keV.在此非平衡態(tài)情況下,光子與電子的相互作用不僅需考慮電子的熱運(yùn)動速度,而且必需在相對論體系下考慮光子與電子的相互作用[2].

利用蒙特卡羅方法模擬高溫等離子體中的輻射輸運(yùn)問題時,需要知道任意能量的光子與任意溫度電子的相互作用(吸收、散射)截面,以及發(fā)生散射后出射光子的能量和角度.本文討論光子與電子的散射截面計(jì)算問題,散射后光子和電子的能譜、角度譜計(jì)算將在后續(xù)的研究中討論.

光子與靜止自由電子發(fā)生Compton散射的微分截面Klein-Nishina公式[3]為式中re為電子經(jīng)典半徑;ε為以電子靜止能量為單位的入射光子能量,即ε=hν/(m0c2),其中hν為入射光子能量,m0為靜止電子質(zhì)量,c為光速(約300 cm/sh,sh是本文的時間單位,1 sh=10?8s);?為入射光子與出射光子夾角(偏轉(zhuǎn)角);?為立體角.

對(1)式關(guān)于立體角積分,得到能量為hν的光子與1個靜止自由電子發(fā)生散射的概率(微觀散射總截面)公式[4]

當(dāng)ε較小時,(2)式的計(jì)算容易產(chǎn)生數(shù)值誤差,且該式中的對數(shù)計(jì)算較費(fèi)時.故通常采用如下的Hastings[5]近似公式代替(2)式計(jì)算微觀散射總截面:

式中η= 1+0.222037ε.這 里 稱(3)式 為KNHastings公式.

在特別設(shè)定的坐標(biāo)系(入射光子沿Z軸飛行、出射光子位于XZ平面)下,光子與速度分布函數(shù)為f(u)的各向同性電子發(fā)生散射的微分截面計(jì)算公式[2]為式中ν和ν′分別為入射、出射光子頻率;ζ=cos?=? ·?′為光子偏轉(zhuǎn)角余弦;u為電子速率;為電子的Lorentz因子;μ=cosα為與電子運(yùn)動方向相關(guān)的第1個積分角度變量,其中α為電子運(yùn)動方向與入射光子方向(Z軸)的夾

角;ε=hν/(m0c2);ε′=hν′/(m0c2);函數(shù)g(u,μ)定義為

?為電子運(yùn)動的方位角.δ函數(shù)是為了滿足散射的動量守恒引入的,且僅當(dāng)

時,(5)式中的δ可以消除掉.電子速度分布f(u)可取任意概率密度函數(shù),這里取相對論麥克斯韋-玻爾茲曼(relativistic Maxwellian-Boltzman,RMB)速度分布函數(shù)(詳見第3節(jié)).

關(guān)于(4)式的計(jì)算比較復(fù)雜,已有多篇文獻(xiàn)介紹如何計(jì)算該式.這些文獻(xiàn)[6?15]的基本思路大多是:首先將(4)式中的積分核函數(shù)展開(冪級數(shù)、勒讓得函數(shù)),然后采用多重?cái)?shù)值積分的方法計(jì)算散射截面以及在此基礎(chǔ)上計(jì)算確定論輸運(yùn)模擬方法[16]所需的分群能譜與角度譜、或者蒙特卡羅輸運(yùn)模擬所需的散射后出射光子的能量與方向[17,18].

本文不討論數(shù)值計(jì)算(4)式,而是利用蒙特卡羅方法,首先逐個直接模擬光子與相對論電子的相互作用,得到若干光子與單個電子散射的概率值,然后取這些概率值的統(tǒng)計(jì)平均值,最后得到一定頻率的光子與一定溫度的電子發(fā)生散射的概率(截面).因此,本文方法的實(shí)質(zhì)是利用蒙特卡羅方法計(jì)算多維復(fù)雜核函數(shù)的積分問題(注:計(jì)算散射截面為6維,計(jì)算散射后的能譜角度譜為8維).

2 光子與相對論電子散射的蒙特卡羅方法

圖1 實(shí)驗(yàn)室系光子與溫度為Te的電子散射示意圖Fig.1.Photon scattering with electrons at temperature Tein lab coordinate.

問題描述:如圖1所示,設(shè)實(shí)驗(yàn)室坐標(biāo)系X0Y0Z0(簡稱0系)中,t時刻有一光子(hν0,?0)(光子能量和飛行方向),求其與r0處溫度為Te的電子發(fā)生散射的概率(微觀散射截面)及散射后光子和電子的能量、角度分布.由于篇幅的原因,本文只討論散射概率的計(jì)算方法,散射后光子和電子的能量、角度分布計(jì)算方法將另文討論.

蒙特卡羅求解步驟:

步驟1根據(jù)電子溫度Te抽樣單個電子狀態(tài)(u0e,?0e)(電子速率和飛行方向).其中,電子飛行方向?0e(?x0e,?y0e,?z0e)為各向同性抽樣:

式中ξ1和ξ2為[0,1]之間的隨機(jī)數(shù).

由此可以求得光子與電子的飛行方向夾角余弦cosα0=?0·?0e.

電子速率u0e從RMB速率分布函數(shù)中抽樣得到,后面將針對此問題單獨(dú)展開討論.

入射光子與抽樣出的單個電子在實(shí)驗(yàn)室系中的散射前狀態(tài)如圖2所示.

圖2 實(shí)驗(yàn)室系中光子與抽樣出的單個電子散射前示意圖Fig.2.Photon scattering with the sampled electron in lab coordinate.

步驟2設(shè)置新笛卡兒坐標(biāo)系X1Y1Z1(簡稱1系),其坐標(biāo)原點(diǎn)與r0重合,X1軸正方向與電子運(yùn)動方向?0e一致,Z1軸垂直于電子運(yùn)動方向與光子飛行方向所組成的平面.由于0系與1系只存在平移和旋轉(zhuǎn)變換,因此在1系中,光子與X1軸(電子飛行方向)夾角α1=α0,光子在1系中的(碰撞前)狀態(tài)為

電子在1系中的(碰撞前)狀態(tài)為

入射光子與單個電子在坐標(biāo)系1中的散射前狀態(tài)如圖3所示.

圖3 坐標(biāo)系1中光子電子散射前示意圖Fig.3.Photon scattering with electron in coordinate 1.

步驟3建立一個隨電子運(yùn)動的笛卡兒坐標(biāo)系X2Y2Z2(簡稱2系).坐標(biāo)系2與坐標(biāo)系1在t時刻重合,運(yùn)動方向?yàn)閄1軸正方向,運(yùn)動速度為u1e(圖4).故電子在2系中始終靜止(u2e=0)且位于坐標(biāo)原點(diǎn).

圖4 坐標(biāo)系1與坐標(biāo)系2關(guān)系示意圖Fig.4.Relationship between coordinate 1 and coordinate 2.

現(xiàn)求光子在2系中的狀態(tài).根據(jù)狹義相對論,對于光子,(Px,Py,Pz,iE/c)構(gòu)成一個四度向量,又因光子動量光子能量E=hν,所以,(ν?x,ν?y,ν?z,iν)亦構(gòu)成四度向量.1系中的四度向量(ν1?x1,ν1?y1,ν1?z1,iν1)與運(yùn)動坐標(biāo)系2系中的四度向量(ν2?x2,ν2?y2,ν2?z2,iν2)可以通過Lorentz變換建立關(guān)系[19]:式中λ=[1?(u1e/c)2]?1/2.

求解(10)式后得到

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

由此計(jì)算得到光子在2系的頻率ν2及飛行方向?2(?x2,?y2,?z2),(12)式中α2是光子在2系中的飛行方向與X2軸的夾角(如圖5所示,Z1,Z2軸均垂直于紙面).

圖5 光子在1系和2系中的飛行方向與坐標(biāo)軸夾角示意圖Fig.5.Photon’s flight directions and angles in coordinate 1 and coordinate 2.

步驟4建立笛卡兒坐標(biāo)系X3Y3Z3(簡稱3系).Z3軸與Z2軸相同,X3軸由X2軸以Z2軸為軸旋轉(zhuǎn)角度α2得到,即X3軸與光子在2系中的飛行方向?2重合,2系與3系的關(guān)系及光子飛行方向如圖6所示,其中的Z軸重合且均垂直于紙面.由于2系與3系之間僅存在旋轉(zhuǎn)變換,故光子頻率不變,因此光子在3系中的狀態(tài)為:hν3=hν2,?x3=1,?y3=0,?z3=0.

至此,在3系中,光子與電子的相互作用圖像是:一個能量為hν3的光子沿X3軸正方向飛行,在坐標(biāo)原點(diǎn)處與一個自由靜止的電子發(fā)生散射,如圖7所示.故可以方便地采用Klein-Nishina公式處理,這里的微觀散射總截面計(jì)算采用(3)式,式中的ε=hν3/(m0c2), 計(jì)算出的σs(ε)就是光子(hν0,?0)與相對論電子(u0e,?0e)發(fā)生散射的概率.

圖6 坐標(biāo)系2與坐標(biāo)系3關(guān)系示意圖Fig.6.Relationship between coordinate 2 and coordinate 3.

圖7 坐標(biāo)系3中光子與電子散射示意圖Fig.7.Photon scattering with electron in coordinate 3.

上述4個步驟描述了入射光子與單個特定(抽樣出的)相對論電子的微觀散射截面計(jì)算方法.由于相對論效應(yīng)引起的空間收縮,與光子飛行方向相反的電子密度比與光子飛行方向相同的電子密度高,因此,需要將截面σs(ε)乘以空間尺度變化引入的Lorentz因子[20].對于光子與溫度為Te的電子發(fā)生散射的概率,重復(fù)步驟1—步驟4若干(N)次,便可以得到N次平均的結(jié)果:

式中σs(ε(n))為第n抽樣計(jì)算得到的微觀散射截面;γ(n)為第n次抽樣計(jì)算的空間變化Lorentz因子:

式中cosα0=?0·?0e為光子與電子的飛行方向夾角余弦.

抽樣足夠多的樣本計(jì)算后,?σs(hν0,Te)?將逐漸收斂到某一固定值,該值即可認(rèn)為是光子(hν0,?0)與溫度為Te的電子發(fā)生散射的概率(微觀散射截面).

分析上面的計(jì)算流程可以看出,對于散射截面,本文方法的計(jì)算精度取決于兩個因素:電子速率分布抽樣的準(zhǔn)確性;抽樣的樣本數(shù)N(收斂性).

3 關(guān)于電子速率的抽樣

經(jīng)典麥克斯韋-玻爾茲曼(classical Maxwellian-Boltzman,CMB)速率分布的概率密度函數(shù)為

此函數(shù)的一種比較新且效率和精度均較高的抽樣方法可以參考文獻(xiàn)[21].

高溫全電離等離子體中,當(dāng)電子達(dá)到熱力學(xué)平衡態(tài)時,溫度為Te的電子速率u服從RMB速率分布,其概率密度函數(shù)[19,22]為

式中k為玻爾茲曼常數(shù)為電子的Lorentz因子;K2(x)為二階變型貝塞爾函數(shù):

求(16)式的積分非常困難,因此抽樣計(jì)算速率u的值很難實(shí)現(xiàn).筆者也未能查閱到關(guān)于(16)式抽樣方法的相關(guān)文獻(xiàn).因此,為了解決此問題,本文推導(dǎo)了一種近似RMB(similar relativistic Maxwellian-Boltzman,SRMB)速率分布函數(shù):

式中β=λ?1.關(guān)于分布函數(shù)的推導(dǎo)及抽樣方法不是本文的重點(diǎn),這里暫且不討論.

本文利用(18)式作為電子速率分布的抽樣函數(shù)分析抽樣結(jié)果.

首先,分析需要多少樣本才能使得抽樣的電子速度分布收斂.針對(18)式編制抽樣程序,計(jì)算并統(tǒng)計(jì)不同抽樣樣本條件下的電子速率分布.圖8給出了在電子溫度為25 keV情況下,樣本數(shù)分別為104,105和106條件下,統(tǒng)計(jì)平均后得到的電子速率概率密度分布及與(18)式解析計(jì)算結(jié)果的比較情況.可以看出當(dāng)抽樣樣本數(shù)達(dá)到106時,抽樣統(tǒng)計(jì)結(jié)果與解析結(jié)果基本一致,由此可以認(rèn)為結(jié)果已收斂.為了盡量排除抽樣樣本引入的誤差,后面的計(jì)算結(jié)果都是樣本數(shù)為107的情況下得到的.

圖8 三種樣本數(shù)條件下的電子速率分布及與解析結(jié)果比較(Te=25 keV)Fig.8.Electron speed distributions and analysis results with three kinds of sampled particles(Te=25 keV).

然后,分析近似(18)式抽樣的電子速率分布的準(zhǔn)確性.這里的準(zhǔn)確性,指的是與RMB速率分布函數(shù)(16)式的符合程度,若符合得越好,則準(zhǔn)確性越高.圖9給出5種電子溫度(1,10,25,50和100 keV)情況下,利用近似(18)式抽樣計(jì)算得到電子速率分布(SRMB sampled),將其與(16)式的解析結(jié)果(RMB analysis)進(jìn)行比較,同時圖中還給出了根據(jù)(15)式抽樣計(jì)算的電子速率分布(CMB sampled).

1)電子溫度低于10 keV情況(圖9(a)和圖9(b))下,相對論效應(yīng)不明顯,利用CMB分布函數(shù)、SRMB分布函數(shù)抽樣計(jì)算得到的結(jié)果與RMB分布函數(shù)的解析結(jié)果基本一致.

2)電子溫度升高至25 keV情況(圖9(c))下,由于相對論效應(yīng),利用CMB分布函數(shù)抽樣計(jì)算得到的結(jié)果與RMB分布函數(shù)的解析結(jié)果有較明顯的差異;然而,利用SRMB分布函數(shù)抽樣計(jì)算得到的結(jié)果與RMB分布函數(shù)的解析結(jié)果仍有較好的符合度.

圖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)隨著電子溫度的繼續(xù)上升 (圖9(d)和圖9(e)),CMB分布函數(shù)抽樣計(jì)算得到的結(jié)果與RMB分布函數(shù)的解析結(jié)果差異繼續(xù)增大且非常顯著;利用SRMB分布函數(shù)抽樣計(jì)算得到的結(jié)果與RMB分布函數(shù)的解析結(jié)果有差異但小于CMB抽樣結(jié)果.另外,從圖9(d)和圖9(e)可以看出,CMB分布抽樣出的電子速度有顯著超光速(c=300 cm/sh)不合理現(xiàn)象出現(xiàn),然而,SRMB分布抽樣出的電子速度不會出現(xiàn)超光速現(xiàn)象.

綜合以上計(jì)算結(jié)果及分析,可得:樣本數(shù)106以上,抽樣計(jì)算結(jié)果基本收斂;電子溫度不超過25 keV,本文的電子速率分布抽樣較準(zhǔn)確,電子溫度高于25 keV,偏差增大,這可能會影響到散射截面的計(jì)算精度.

4 散射截面的數(shù)值模擬結(jié)果與分析

根據(jù)本文第2節(jié)介紹的方法,研制了光子-相對論麥克斯韋電子散射模擬程序,該程序可以計(jì)算任意光子(能量范圍為0—1 MeV)與任意溫度電子(溫度范圍為0—100 keV)的微觀散射截面及散射角度譜、能譜.表1是文獻(xiàn)[11]中給出的利用數(shù)值積分方法計(jì)算的光子-電子散射截面.為了分析電子熱運(yùn)動對散射截面影響,表2列出了利用(3)式計(jì)算的光子與靜止自由電子的散射截面;表3是利用本文方法計(jì)算的光子-電子散射截面.為了便于比較,表2和表3在給出散射截面的同時還列出了其與表1中對應(yīng)值的相對誤差其中σtable1,σtable2和σtable3分別是表1、表2和

首先,分析表1可以看出:光子能量越高,截面越小;電子熱運(yùn)動使得截面變小,且電子溫度越高、熱運(yùn)動平均速度越大,截面減小越多.

表1 光子-相對論麥克斯韋電子散射截面(文獻(xiàn)[11]數(shù)值積分計(jì)算)Table 1.Photon-Maxwellian electron cross-sections(numerical integration method,Ref.[11]).

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

表3 光子-相對論電子散射截面(本文蒙特卡羅方法計(jì)算)Table 3.Photon-Maxwellian electron cross-sections(proposed Monte Carlo method).

表2的截面計(jì)算因未考慮電子熱運(yùn)動,故相對于表1其值偏大,且電子溫度越高,截面偏大越多;光子能量越高,截面偏大越多;例如:在光子能量為1 MeV、電子溫度為100 keV情況下,相對誤差達(dá)18%左右.

表3本文方法計(jì)算的散射截面,在電子溫度不高于25 keV情況下與表1的結(jié)果很接近,相對誤差不超過2%;但是隨著電子溫度繼續(xù)升高,偏差逐漸增大,例如:在光子能量為1 MeV、電子溫度為100 keV情況下,相對誤差也達(dá)到7.8%.

分析本文方法在電子溫度較高(>25 keV)情況下與國外的數(shù)值積分方法差異的原因:前面已經(jīng)提到過,電子平均速度越大截面越小.然而,從前面的圖9(d)和圖9(e)可知,在電子溫度較高時,本文目前關(guān)于電子速率的抽樣方法抽出的電子速率整體偏小,這可能就是導(dǎo)致本文方法計(jì)算的散射截面偏大的原因.因此,本研究下一步將圍繞相對論電子速率的準(zhǔn)確抽樣目標(biāo)開展.

5 結(jié) 論

本文提出了一種計(jì)算光子與相對論麥克斯韋分布電子散射截面的方法.該方法的實(shí)質(zhì)是利用蒙特卡羅方法計(jì)算多維復(fù)雜核函數(shù)的積分.本文方法具有如下特點(diǎn):第一,在給定的光子能量范圍(0—1 MeV)和電子溫度范圍(0—100 keV)內(nèi),不需要針對光子能量和電子溫度分群離散,即可以根據(jù)單一光子能量和單一電子溫度直接計(jì)算出對應(yīng)的散射截面;第二,散射后光子與電子的能量和角度是連續(xù)的,可以方便地利用蒙特卡羅方法確定(另文介紹),因此,本文方法特別適用于利用蒙特卡羅方法模擬光子在高溫等離子體中的Compton散射和逆Compton散射問題;第三,與數(shù)值積分方法計(jì)算散射截面相比,本文方法無需做復(fù)雜繁瑣的數(shù)學(xué)處理,因此更加直觀、簡潔.

本文方法在電子溫度不超過25 keV情況下的計(jì)算精度較高(注:此情況能夠滿足大多數(shù)ICF、熱核裝置等問題所需);在電子溫度更高情況下,可能是目前的電子速率抽樣有偏差的原因,導(dǎo)致目前的計(jì)算結(jié)果與國外多重?cái)?shù)值積分結(jié)果存在較明顯差異;待將來找到更好的相對論電子速率抽樣方法后,此缺陷應(yīng)該能夠被克服.