基于粒子云方法研究地磁場(chǎng)對(duì)運(yùn)動(dòng)電子束的時(shí)空分布影響

馮明航，付松，倪彬彬，雷良建，焦鹿懷

（武漢大學(xué)電子信息學(xué)院，武漢 430072）

0 引言

地球空間存在大量的帶電粒子束結(jié)構(gòu)[1-2]，一般而言，自然的電子束往往可以由波粒相互作用產(chǎn)生的粒子沉降或者等離子體不穩(wěn)定性引起[3-6]。這些粒子束中的帶電粒子往往具有很高的動(dòng)能和很強(qiáng)的方向性，因此不但可以影響空間天氣和空間環(huán)境[7]，還可能通過(guò)總劑量輻射損傷、單粒子效應(yīng)和充放電效應(yīng)等對(duì)運(yùn)行在地球空間的人造天體造成潛在的重要危害[8-10]。長(zhǎng)期以來(lái)，空間環(huán)境中的帶電粒子束尤其是電子束一直是空間物理學(xué)的重要研究?jī)?nèi)容，科學(xué)家通過(guò)天基、地基觀測(cè)，理論分析以及數(shù)值模擬來(lái)研究電子束在空間環(huán)境中的形成、演化以及消亡等關(guān)鍵過(guò)程[11-12]，希望能夠準(zhǔn)確預(yù)測(cè)電子束的時(shí)空分布來(lái)降低或者避免其對(duì)人造天體的危害。

研究電子束在地球空間環(huán)境所起作用的一個(gè)重要問(wèn)題是，電子束在傳輸過(guò)程中自場(chǎng)效應(yīng)以及對(duì)背景地磁場(chǎng)的響應(yīng)的空間電荷效應(yīng)耦合結(jié)果，將如何影響電子束的結(jié)構(gòu)形態(tài)以及電子束本身的相空間演化[13]。一般而言，電子束在運(yùn)動(dòng)過(guò)程中，帶電粒子相互之間受庫(kù)侖斥力影響，會(huì)產(chǎn)生橫向結(jié)構(gòu)上的發(fā)散以及縱向結(jié)構(gòu)上的拉伸；此外帶電粒子受外在電磁場(chǎng)所產(chǎn)生的洛倫茲力影響，會(huì)沿外在電場(chǎng)方向加速并圍繞外在磁場(chǎng)發(fā)生偏轉(zhuǎn)[14]。以上多種效應(yīng)耦合在一起，構(gòu)成了研究電子束空間電荷效應(yīng)的物理機(jī)制背景。針對(duì)此復(fù)雜的多物理機(jī)制耦合問(wèn)題，常用的研究方法有：1）基于電荷中和的電子束傳輸線性解析理論，以求解電子束包絡(luò)方程為手段研究電子束空間電荷效應(yīng)的演化問(wèn)題[15]；2）基于單粒子運(yùn)動(dòng)方程的試驗(yàn)粒子方法，以求解帶電粒子的單粒子運(yùn)動(dòng)方程來(lái)描述典型粒子的運(yùn)動(dòng)軌跡[16]。方法1 對(duì)于電子束傳輸?shù)牡碗A線性和單一物理過(guò)程可以給出非常直觀的物理解釋?zhuān)鐭o(wú)背景場(chǎng)的電子束空間電荷效應(yīng)問(wèn)題，但是對(duì)于研究非線性物理過(guò)程和多物理機(jī)制耦合過(guò)程則比較乏力；方法2可以很好地描述單一參考粒子對(duì)外加電磁場(chǎng)的運(yùn)動(dòng)響應(yīng)，但是無(wú)法體現(xiàn)電子束內(nèi)部帶電粒子之間相互作用的效果。針對(duì)以上兩種方法的優(yōu)勢(shì)和不足，需要使用Vlasov-Poisson 動(dòng)理論來(lái)研究電子束的空間電荷效應(yīng)問(wèn)題[17]；而粒子云（particle-in-cell,PIC）方法是目前最為有效的解決該問(wèn)題的數(shù)值模擬解法[18-19]，能夠完整自洽地描述Vlasov 系統(tǒng)中的粒子相空間以及電磁場(chǎng)時(shí)空分布變化，并被廣泛應(yīng)用于空間物理和天體物理、受控?zé)岷司圩?、等離子體物理和真空電子學(xué)等多個(gè)領(lǐng)域[20-21]。

為了研究電子束空間電荷效應(yīng)在地球空間環(huán)境中的作用及影響，我們開(kāi)發(fā)了靜電模式的PIC代碼，并建立了基于PIC算法的數(shù)值模型。本文給出了該模型的核心算法、電磁場(chǎng)求解方法、帶電粒子運(yùn)動(dòng)方程求解方法和數(shù)值算法初始化等，并對(duì)該算法的可靠性進(jìn)行驗(yàn)證；進(jìn)而采用該模型對(duì)不同地磁條件下不同能量的電子束的空間電荷效應(yīng)進(jìn)行定量分析及對(duì)比。

1 理論模型及數(shù)值方法

為了定量研究地球空間環(huán)境帶電粒子束受自場(chǎng)和背景地磁場(chǎng)的聯(lián)合效應(yīng)，需要同時(shí)考慮由庫(kù)侖力產(chǎn)生的帶電粒子自場(chǎng)效應(yīng)和由洛倫茲力產(chǎn)生的帶電粒子對(duì)背景地磁場(chǎng)的響應(yīng)，因此使用描述帶電粒子動(dòng)力學(xué)過(guò)程的自洽PIC方法來(lái)進(jìn)行數(shù)值建模。

1.1 粒子云模型算法

PIC算法被廣泛用來(lái)研究帶電粒子束空間電荷效應(yīng)產(chǎn)生的集體效應(yīng)，該方法能同時(shí)處理帶電粒子分布對(duì)電磁場(chǎng)的影響以及電磁場(chǎng)對(duì)帶電粒子運(yùn)動(dòng)的影響。為此我們開(kāi)發(fā)了針對(duì)帶電粒子束運(yùn)動(dòng)的靜電模式PIC 代碼（以下簡(jiǎn)稱(chēng)beamPIC）。如圖1所示，算法的主循環(huán)包含4個(gè)步驟[22]：1）使用CIC（cloudin-cell）方法將電荷分配至離散化的空間網(wǎng)格格點(diǎn)，再求得各空間格點(diǎn)上的電荷密度分布；2）根據(jù)網(wǎng)格格點(diǎn)處的電荷密度分布，使用基于快速傅里葉變換（FFT）的積分形式格林函數(shù)法求解Poisson 方程，從而得到空間網(wǎng)格點(diǎn)上的電勢(shì)分布，再對(duì)電勢(shì)進(jìn)行差分得到網(wǎng)格點(diǎn)上的電場(chǎng)分布；3）利用插值獲得每個(gè)粒子位置處的電磁場(chǎng)；4）求解洛倫茲方程，利用Boris方法更新粒子速度，再使用蛙跳法求得粒子的坐標(biāo)。

圖1 粒子云模型的核心算法Fig.1 Core algorithm of PIC method

1.2 空間電磁場(chǎng)求解

常用的地球磁場(chǎng)模型包括：

1）以國(guó)際地磁學(xué)和高層大氣物理學(xué)協(xié)會(huì)（IAGA）頒布的國(guó)際地磁參考場(chǎng)模型（IGRF）。該模型使用球諧函數(shù)擬合地球的內(nèi)源場(chǎng)；每5年頒布一個(gè)新的模型，目前最新版本是2019年12月發(fā)布的IGRF-13模型[23]。

2）以Tsyganenko模型為代表的半經(jīng)驗(yàn)?zāi)Ｐ?。該模型根?jù)衛(wèi)星的磁場(chǎng)觀測(cè)數(shù)據(jù)和物理機(jī)制建立，將地球磁場(chǎng)看作內(nèi)源場(chǎng)和外源場(chǎng)之和。該模型從發(fā)布之后不斷地添加新的衛(wèi)星數(shù)據(jù)來(lái)更新模型參數(shù)，目前較為廣泛使用的版本包括T96和TS05等[24]。

為了簡(jiǎn)化問(wèn)題，本研究所使用的空間電磁場(chǎng)包括了由三維偶極子場(chǎng)模型[25]得到的背景地磁場(chǎng)和帶電粒子束的空間電荷效應(yīng)引起的自場(chǎng)兩部分。背景地磁場(chǎng)為

式中：?為電勢(shì)；G為格林函數(shù)；ρ為電荷密度分布函數(shù)；(x,y,z)為場(chǎng)點(diǎn)坐標(biāo)；(x',y',z')為源點(diǎn)坐標(biāo)。由于本文研究自由空間的帶電粒子束在地磁場(chǎng)下的動(dòng)力學(xué)過(guò)程，所以選取開(kāi)放邊界條件，此時(shí)G為

設(shè)定模擬區(qū)域?yàn)槿S自由空間中的長(zhǎng)方體，其在x、y、z方向上的長(zhǎng)度分別為L(zhǎng)x、Ly、Lz，同時(shí)在三個(gè)方向上的節(jié)點(diǎn)數(shù)分別為Nx、Ny、Nz，在這些節(jié)點(diǎn)上放置點(diǎn)電荷，則任意網(wǎng)格點(diǎn)上的電勢(shì)為

式中：hx、hy和hz分別為x、y和z方向上的網(wǎng)格尺寸；xi=(i?1)hx，yj=(j?1)hy，zk=(k?1)hz；對(duì)于奇異值點(diǎn)（i=i'，j=j'，k=k'），G取值為1。利用上式求得?(x,y,z)后，再通過(guò)差分求解E=-??即可得到該節(jié)點(diǎn)上的電場(chǎng)強(qiáng)度。對(duì)于式(4)，直接利用卷積數(shù)值求和，時(shí)間復(fù)雜度為N6（N為網(wǎng)格點(diǎn)總數(shù)）。為了利用FFT算法本身的周期性特點(diǎn)，并更好地滿(mǎn)足開(kāi)放邊界條件，可以使用Hockney Trick 方法[19]，將計(jì)算域擴(kuò)展為2Nx×2Ny×2Nz，同時(shí)將格林函數(shù)和電荷密度分布函數(shù)也拓展至新的計(jì)算域空間，計(jì)算得到電勢(shì)后，取實(shí)際物理場(chǎng)域Nx×Ny×Nz的值即為空間電荷產(chǎn)生的電勢(shì)。該方法可將時(shí)間復(fù)雜度優(yōu)化至N4lgN。

由上可知，通過(guò)使用三維積分形式的格林函數(shù)法，可以得到非常準(zhǔn)確的空間電場(chǎng)分布，具體實(shí)現(xiàn)步驟為：1）根據(jù)空間分布范圍確定網(wǎng)格尺寸參數(shù)；2）利用三維線性加權(quán)CIC 方法將電荷分配到鄰近的8個(gè)網(wǎng)格點(diǎn)上，從而得到 ρD；3）利用三維積分格林函數(shù)解析方法求得GD；4)對(duì) ρD和GD作三維FFT 變換，得到和，進(jìn)而求得；5）對(duì) ??D作三維FFT 逆變換，得到網(wǎng)格點(diǎn)上的 ?D；6）利用數(shù)值差分求解ED=-??D得到空間電場(chǎng)分布；7）利用插值方法將網(wǎng)格點(diǎn)上的電場(chǎng)插值到粒子所在位置，得到粒子位置處的電場(chǎng)強(qiáng)度。

通過(guò)以上數(shù)值方法求解Poisson 方程，可以求得電子束空間電荷效應(yīng)產(chǎn)生的電場(chǎng)分布。

1.3 粒子運(yùn)動(dòng)方程求解

由上述方法得到背景地磁場(chǎng)和電子束的自場(chǎng)空間分布，通過(guò)帶電粒子受力分析可以求得電子穿越電磁場(chǎng)時(shí)所受的洛倫茲力。采用數(shù)值方法求解描述電子運(yùn)動(dòng)的相對(duì)論洛倫茲方程，使用具有二階精度的蛙跳格式來(lái)推進(jìn)粒子運(yùn)動(dòng)：首先，粒子運(yùn)動(dòng)半個(gè)時(shí)間步，再通過(guò)求解電子束靜止坐標(biāo)系下的泊松方程計(jì)算空間電荷效應(yīng)；通過(guò)洛倫茲變換獲得實(shí)驗(yàn)室坐標(biāo)系下的電磁場(chǎng)，然后采用Boris方法[26]更新粒子動(dòng)量；之后粒子再運(yùn)動(dòng)半個(gè)時(shí)間步，重復(fù)上述步驟，直至滿(mǎn)足停止條件。

帶電粒子的運(yùn)動(dòng)學(xué)方程形式如下：

其中Δt為設(shè)置的時(shí)間步長(zhǎng)。Boris方法廣泛應(yīng)用于數(shù)值等離子體模擬中，其具有二階精度和簡(jiǎn)潔的數(shù)值形式，式(7)和式(10)計(jì)算了電場(chǎng)力對(duì)粒子動(dòng)量的作用，式(8)和式(9)計(jì)算了磁場(chǎng)對(duì)粒子動(dòng)量的改變。

1.4 算法初始化與終止條件

使用PIC程序，通過(guò)讀取配置文件進(jìn)行初始化，參數(shù)包含模擬系統(tǒng)環(huán)境設(shè)置和電子束生成參數(shù)，其中：系統(tǒng)環(huán)境函數(shù)可以設(shè)定模擬的時(shí)間步長(zhǎng)、最大時(shí)間步數(shù)、模擬空間的網(wǎng)格數(shù)、電子束運(yùn)動(dòng)距離以及背景磁場(chǎng)等參數(shù)；電子束生成函數(shù)包含粒子束分布類(lèi)型、能量、宏粒子數(shù)量、電荷量以及質(zhì)量等參數(shù)。對(duì)于相空間的電子束，使用電子束尺寸和發(fā)射度兩個(gè)概念來(lái)描述電子束在生成和運(yùn)動(dòng)過(guò)程中的相空間演化情況，包括均方根（RMS）尺寸和均方根偏轉(zhuǎn)角[27]。對(duì)于電子束的初始分布，目前主要研究Uniform 分布和Gauss分布等情況。

程序在完成初始化和檢查配置文件無(wú)誤后開(kāi)始進(jìn)入PIC算法的循環(huán)計(jì)算。當(dāng)電子束滿(mǎn)足設(shè)定的終止條件如運(yùn)動(dòng)最遠(yuǎn)距離或時(shí)間步數(shù)達(dá)到最大值時(shí)程序終止循環(huán)。程序運(yùn)行的每個(gè)時(shí)間上都可將粒子的位置空間和速度空間的信息記錄下來(lái)得到完整的電子束演化過(guò)程。

2 模擬結(jié)果及分析

2.1 初始設(shè)定

圖2所示為Gauss分布和Uniform 分布的兩類(lèi)電子束在初始時(shí)刻的橫向尺寸（圖上）和偏轉(zhuǎn)角（圖下）分布，其中：橫向均方根尺寸為5 mm，橫向均方根偏轉(zhuǎn)角為0.001 mrad，縱向尺寸為5 mm，縱向能散度為1%，每次模擬所使用的宏粒子總數(shù)均為50萬(wàn)個(gè)。在后續(xù)模擬過(guò)程中均使用上述參數(shù)生成電子束。

圖2 電子束Gauss 初始分布和Uniform 初始分布的橫向尺寸和偏轉(zhuǎn)角分布Fig.2 Transverse RMS size and rotation angle of electron beamsfor initial Gaussand Uniform distributions

2.2 可靠性分析

為了進(jìn)一步驗(yàn)證所開(kāi)發(fā)PIC算法的可靠性，將beamPIC的模擬結(jié)果與ASTRA（A Space Charge Tracking Algorithm）軟件包的進(jìn)行對(duì)比分析。ASTRA[28]是由德國(guó)電子同步加速器研究所發(fā)布的空間電荷追蹤算法程序，包含了完整的粒子初始分布生成算法和粒子追蹤算法等，能夠通過(guò)PIC算法來(lái)計(jì)算帶電粒子的空間電荷效應(yīng)等多種集體自場(chǎng)效應(yīng)，被廣泛應(yīng)用于高能帶電粒子的仿真研究。需要說(shuō)明的是，ASTRA 的算法中假設(shè)電子束在縱向具有無(wú)限長(zhǎng)度，因此只能計(jì)算該類(lèi)電子束的橫向四維相空間變化，而我們的beamPIC算法能夠研究有限長(zhǎng)度電子束的六維相空間變化。

為了對(duì)兩種方法進(jìn)行對(duì)比，將能量為100 keV的電子束同時(shí)輸入到beamPIC程序和ASTRA 軟件，在相同環(huán)境參數(shù)設(shè)定下，比較電子束傳輸100 m過(guò)程中均方根橫向尺寸和均方根偏轉(zhuǎn)角變化情況，結(jié)果如圖3所示，其中：藍(lán)色實(shí)線和紅色實(shí)線分別為初始Gauss分布電子束的beamPIC和ASTRA模擬結(jié)果，黃色虛線和紫色虛線分別為初始Uniform 分布電子束的beamPIC和ASTRA 模擬結(jié)果。對(duì)于背景磁場(chǎng)，選取了無(wú)背景磁場(chǎng)（圖3(a)）、與電子束發(fā)射方向平行的背景磁場(chǎng)（圖3(b)）以及與電子束發(fā)射方向垂直的背景磁場(chǎng)（圖3(c)）3種情況，地磁場(chǎng)大小取L-shell為5時(shí)的磁赤道面磁場(chǎng)大?。ㄓ墒?1)可知為249.6 nT）。

圖3 電子束傳輸?shù)腷eamPIC與ASTRA 結(jié)果對(duì)比Fig.3 Transfer of electron beams obtained by using beamPIC and ASTRA methods

從圖3可以看出：當(dāng)電子束運(yùn)動(dòng)距離為100 m時(shí)，針對(duì)所研究的3種背景場(chǎng)和2種初始分布設(shè)定，beamPIC程序和ASTRA 程序結(jié)果一致性都非常好，RMS橫向尺寸和RMS偏轉(zhuǎn)角變化基本完全吻合；電子束運(yùn)動(dòng)至100 m 時(shí)，RMS橫向尺寸約為230 mm，RMS偏轉(zhuǎn)角約為2.3 mrad。上述結(jié)果說(shuō)明對(duì)比ASTRA 軟件，beamPIC程序計(jì)算結(jié)果具有較好的一致性，對(duì)于各種初始分布和各種背景磁場(chǎng)情況均能保持很好的數(shù)值精度，具有良好的可靠性。

2.3 參量化分析

為了進(jìn)一步研究不同地磁場(chǎng)影響下電子束運(yùn)動(dòng)的空間電荷效應(yīng)，利用beamPIC程序計(jì)算初始Gauss分布和Uniform 分布的兩類(lèi)電子束在不同大小背景磁場(chǎng)作用下運(yùn)動(dòng)100 m 時(shí)的RMS橫向尺寸、偏轉(zhuǎn)角和縱向尺寸變化，來(lái)分析電子束在地球磁場(chǎng)和空間電荷效應(yīng)耦合下的相空間變化情況。為了方便計(jì)算，此處只研究了電子束縱向垂直于背景磁場(chǎng)的情況。

在beamPIC的參量計(jì)算中，L-shell 分別選取3、4、5、6和7，根據(jù)Dipole場(chǎng)模型，計(jì)算出背景磁場(chǎng)大小分別為1155.6 nT、487.5 nT、249.6 nT、144.4 nT和91.0 nT。電子束能量分別選取0.1 keV、1 keV、10 keV 和100 keV，圖4展示了各能量電子在垂直于背景磁場(chǎng)時(shí)的回旋半徑?？梢园l(fā)現(xiàn)0.1 keV 電子在L-shell 為3和4，以及1 keV 電子在L-shell 為3時(shí)回旋半徑小于傳輸距離100 m，因此在下面的計(jì)算中忽略相應(yīng)的這3種情況。

圖4 0.1、1、10和100 keV 電子在地磁場(chǎng)下的回旋半徑Fig.4 Cyclotron radius of electrons of energies of 0.1,1,10 and 100 keV in the geomagnetic field

2.3.1 L-shell=3的模擬結(jié)果

圖5所示為L(zhǎng)-shell=3時(shí)10 keV 和100 keV 的電子束運(yùn)動(dòng)100 m 的RMS橫向尺寸、偏轉(zhuǎn)角和RMS縱向尺寸變化情況。其中：黃線表示10 keV 電子束的結(jié)果，紫線表示100 keV 電子束的結(jié)果；實(shí)線表示Gauss初始分布電子束，虛線表示Uniform 初始分布電子束?？梢钥闯鰞煞N初始分布的結(jié)果在100 m運(yùn)動(dòng)距離上比較接近。再以10 keV 電子束為例，從圖5可以看出運(yùn)動(dòng)100 m 距離時(shí)，RMS橫向尺寸約為1 m，RMS偏轉(zhuǎn)角約為10 mrad，RMS縱向尺寸約為1.2 m。其中，RMS橫向和縱向尺寸的變化主要由電子束自場(chǎng)效應(yīng)所產(chǎn)生的擴(kuò)散和拉伸引起，而RMS偏轉(zhuǎn)角的變化主要由對(duì)背景場(chǎng)的偏轉(zhuǎn)響應(yīng)和能量色散引起。此外還可以看到，100 keV 電子束橫向和縱向發(fā)散明顯小于10 keV 電子束，這說(shuō)明電子束能量越高，自場(chǎng)產(chǎn)生的發(fā)散和拉伸效應(yīng)以及受背景場(chǎng)的影響所產(chǎn)生的電子束相空間的變化越小。

圖5 L-shell=3時(shí)的模擬結(jié)果（B0=1155.6 nT）Fig.5 Simulation results for the case of L-shell = 3(B0 =1155.6 nT)

2.3.2 L-shell=4的模擬結(jié)果

圖6所示為L(zhǎng)-shell=4時(shí)的1 keV、10 keV 和100 keV 電子束運(yùn)動(dòng)100 m 的RMS橫向尺寸、偏轉(zhuǎn)角和縱向尺寸變化情況，其中紅色實(shí)線和虛線分別表示初始Gauss分布和初始Uniform 分布的1 keV電子束的結(jié)果。運(yùn)動(dòng)100 m 距離時(shí)，1 keV 電子束的RMS橫向尺寸約為4 m，RMS偏轉(zhuǎn)角約為35 mrad，RMS縱向尺寸約為3.5 m，其橫向、縱向發(fā)散程度和偏轉(zhuǎn)角均遠(yuǎn)大于10 keV 和100 keV 電子束的結(jié)果。此外，隨著背景磁場(chǎng)的變?nèi)酰?0 keV 電子束和100 keV 電子束在橫向和縱向發(fā)散上相較于L-shell為3時(shí)有所減小。總之，對(duì)于較低能量電子束而言，其回旋半徑較小，在運(yùn)動(dòng)相同距離上受背景磁場(chǎng)影響更大，同時(shí)所需時(shí)間更長(zhǎng)，因此橫向和縱向的相空間變化均最為明顯。

圖6 L-shell=4時(shí)模擬結(jié)果（B0=487.5 nT）Fig.6 Simulation results for the case of L-shell = 4(B0 =487.5 nT)

2.3.3 L-shell=5的模擬結(jié)果

圖7所示為L(zhǎng)-shell=5時(shí)的0.1 keV、1 keV、10 keV 和100 keV 電子束運(yùn)動(dòng)100 m 的RMS橫向尺寸、偏轉(zhuǎn)角和縱向尺寸變化情況，其中藍(lán)色實(shí)線和虛線分別表示初始Gauss分布和初始Uniform分布的0.1 keV 電子束的結(jié)果。從圖7可以看到，0.1 keV 電子束運(yùn)動(dòng)100 m 距離時(shí)，RMS橫向尺寸約為17 m，RMS偏轉(zhuǎn)角約為150 mrad，RMS縱向尺寸約為11 m，其橫向、縱向發(fā)散程度和偏轉(zhuǎn)角均遠(yuǎn)大于其他3組能量電子束的結(jié)果。

圖7 L-shell=5時(shí)模擬結(jié)果（B0 = 249.6 nT）Fig.7 Simulation resultsfor the caseof L-shell = 5(B0 =249.6 nT)

2.3.4 L-shell=6的模擬結(jié)果

圖8所示為L(zhǎng)-shell=6時(shí)的0.1 keV、1 keV、10 keV 和100 keV 電子束運(yùn)動(dòng)100 m 的RMS橫向尺寸、偏轉(zhuǎn)角和縱向尺寸變化情況。可以看出：0.1 keV 電子束運(yùn)動(dòng)100 m 距離時(shí)，RMS橫向尺寸約為19 m，RMS偏轉(zhuǎn)角約為112 mrad，RMS縱向尺寸約為10 m，其橫向、縱向發(fā)散程度和偏轉(zhuǎn)角均遠(yuǎn)大于其他3組能量電子束；隨著L-shell的增大，背景磁場(chǎng)減小，4組能量電子束在橫向發(fā)散上均有所減小，在縱向發(fā)散上0.1 keV 電子束有所減小，而其他3組變化不明顯；對(duì)于0.1 keV 電子束，相比于L-shell=5時(shí)的結(jié)果，L-shell=6時(shí)其橫向和縱向發(fā)散均有明顯的減小。

圖8 L-shell等于6時(shí)模擬結(jié)果Fig.8 Simulation resultsfor the case of L-shell = 6(B0=144.4 nT)

2.3.5 L-shell=7的模擬結(jié)果

圖9所示為L(zhǎng)-shell=7時(shí)的0.1 keV、1 keV、10 keV 和100 keV 電子束運(yùn)動(dòng)100 m 的RMS橫向尺寸、偏轉(zhuǎn)角和縱向尺寸變化情況?？梢钥吹剑?.1keV電子束運(yùn)動(dòng)100 m 距離時(shí)，RMS橫向尺寸約為11 m，RMS偏轉(zhuǎn)角約為105 mrad，RMS縱向尺寸約為10 m，其橫向、縱向發(fā)散程度和偏轉(zhuǎn)角均遠(yuǎn)大于其他3組能量電子束；相比于L-shell=6時(shí)的結(jié)果，L-shell=7時(shí)0.1 keV 電子束RMS橫向尺寸和偏轉(zhuǎn)角進(jìn)一步減小，但RMS縱向尺寸幾乎不變，其他3組在橫向和縱向發(fā)散趨勢(shì)上基本保持不變。

圖9 L-shell=7時(shí)的模擬結(jié)果（B0 =91.0 nT）Fig.9 Simulation resultsfor the case of L-shell =7(B0= 91.0 nT)

2.4 RMS橫向尺寸、偏轉(zhuǎn)角及縱向尺寸隨地磁場(chǎng)強(qiáng)度變化結(jié)果對(duì)比

圖10展示了初始為Gauss和Uniform 兩種分布，能量為0.1、1、10和100 keV 的電子束在不同背景磁場(chǎng)作用下運(yùn)動(dòng)100 m 時(shí)的RMS橫向尺寸。0.1 keV 電子束在L-shell=5時(shí)RMS橫向尺寸有最大值約為17 m，在運(yùn)動(dòng)60～100 m 階段，由于背景地磁場(chǎng)影響，橫向尺寸增長(zhǎng)速度變快；1 keV 電子束在L-shell=4時(shí)RMS橫向尺寸取得最大值約3.7 m，其他背景磁場(chǎng)下橫向尺寸最終集中在3.2 m 左右；10 keV 電子束在L-shell=3時(shí)RMS橫向尺寸有最大值約為1 m，最小值在L-shell=7時(shí)，約為0.95 m；100 keV 電子束RMS橫向尺寸在L-shell 從3～7變化過(guò)程中均增長(zhǎng)至約0.23 m，表明其受背景磁場(chǎng)影響較小。

圖11所示為電子束運(yùn)動(dòng)100 m 過(guò)程中RMS偏轉(zhuǎn)角的變化。能量為0.1 keV 電子束在L-shell=5時(shí)偏轉(zhuǎn)最強(qiáng)，最終RMS偏轉(zhuǎn)角約為150 mrad，Lshell為6和7時(shí)的RMS偏轉(zhuǎn)角也快速增長(zhǎng)超過(guò)100 mrad；能量為1 keV 電子束RMS偏轉(zhuǎn)角最大值降低到約36 mrad，此時(shí)L-shell=4，其他L-shell下的RMS偏轉(zhuǎn)角增長(zhǎng)至31 mrad 左右；10 keV 電子束在L-shell=3時(shí)RMS偏轉(zhuǎn)角增長(zhǎng)至最大值約為10.5 mrad，表明其受背景磁場(chǎng)影響進(jìn)一步減小，其他L-shell 下RMS偏轉(zhuǎn)角集中在9.5 mrad 附近；100 keV 電子束在L-shell 變化過(guò)程中RMS偏轉(zhuǎn)角最終均增長(zhǎng)至2.4 mrad 附近。

圖10 電子束在地磁場(chǎng)運(yùn)動(dòng)過(guò)程中RMS橫向尺寸變化結(jié)果Fig.10 Variations of transverse size of electron beam travelling through the geomagnetic field

圖11 電子束在地磁場(chǎng)運(yùn)動(dòng)過(guò)程中RMS偏轉(zhuǎn)角變化結(jié)果Fig.11 Variations of rotation angle of electron beam travelling through the geomagnetic field

圖12所示為電子束運(yùn)動(dòng)100 m 過(guò)程中RMS縱向尺寸的變化。0.1 keV 電子束在L-shell=5時(shí)縱向尺寸有最大值約11 m，其他磁場(chǎng)條件下集中在10 m 左右；1 keV、10 keV 和100 keV 電子束縱向尺寸受背景磁場(chǎng)影響較小，運(yùn)動(dòng)結(jié)束時(shí)分別集中在約3.1 m、1.1 m 和0.5 m，此時(shí)電子束縱向尺寸變化主要受電子束空間電荷效應(yīng)影響。

圖12 電子束在地磁場(chǎng)運(yùn)動(dòng)過(guò)程中RMS縱向尺寸變化結(jié)果Fig.12 Variations of longitudinal size of electron beam travelling through the geomagnetic field

3 結(jié)束語(yǔ)

本文基于PIC方法建立了三維電子束的靜電動(dòng)力學(xué)模擬系統(tǒng)beamPIC程序，并研究了在地球磁場(chǎng)作用下電子束運(yùn)動(dòng)過(guò)程中由于自場(chǎng)和背景場(chǎng)的聯(lián)合效應(yīng)而出現(xiàn)的發(fā)散、拉伸和偏轉(zhuǎn)現(xiàn)象，得到如下結(jié)論：

1）基于PIC方法的電子束自洽模型可以很好地解決帶電粒子與場(chǎng)的耦合過(guò)程，beamPIC程序整體具有良好的可靠性。

2）電子束在空間運(yùn)動(dòng)傳輸100 m 過(guò)程中，受地球磁場(chǎng)和自身空間電荷效應(yīng)耦合作用，會(huì)出現(xiàn)發(fā)散和偏轉(zhuǎn)現(xiàn)象：背景磁場(chǎng)強(qiáng)度一定時(shí)，電子束的RMS橫向尺寸、偏轉(zhuǎn)角和縱向尺寸與電子束能量成負(fù)相關(guān)；電子束能量一定時(shí)，其RMS橫向尺寸、偏轉(zhuǎn)角與背景磁場(chǎng)強(qiáng)度成正相關(guān)；空間電子束能量＞1 keV 時(shí)，其RMS縱向尺寸在不同背景磁場(chǎng)條件下變化基本一致，此時(shí)電子束縱向主要由于空間電荷效應(yīng)而產(chǎn)生發(fā)散。

上述研究結(jié)果表明，電子束結(jié)構(gòu)在背景地磁場(chǎng)環(huán)境中的運(yùn)動(dòng)演化非常復(fù)雜，會(huì)同時(shí)受到束結(jié)構(gòu)自身的靜電力和背景磁場(chǎng)的洛倫茲力作用，使電子束產(chǎn)生橫向的發(fā)散、縱向的拉伸以及繞磁力線的偏轉(zhuǎn)和能量色散等多種效應(yīng)，且各種效應(yīng)之間往往存在復(fù)雜的耦合效應(yīng)。使用粒子云算法可以對(duì)這一復(fù)雜物理過(guò)程進(jìn)行精確有效的描述。

值得注意的是，在本研究中，我們假設(shè)背景空間為真空環(huán)境，而實(shí)際的空間環(huán)境中存在著背景等離子體和各種波動(dòng)現(xiàn)象。因此在未來(lái)的研究中，我們將著重考察電子束結(jié)構(gòu)在運(yùn)動(dòng)過(guò)程中與背景等離子體和波動(dòng)的相互作用。