行星際日冕物質(zhì)拋射引起福布斯下降的一維隨機(jī)微分模擬

倪素蘭顧斌韓智伊

1)(南京信息工程大學(xué)物理系,南京 210044)

2)(南京信息工程大學(xué)空間天氣研究所,南京 210044)

行星際日冕物質(zhì)拋射引起福布斯下降的一維隨機(jī)微分模擬

倪素蘭1)2)顧斌1)2)?韓智伊1)2)

1)(南京信息工程大學(xué)物理系,南京 210044)

2)(南京信息工程大學(xué)空間天氣研究所,南京 210044)

(2017年3月6日收到;2017年3月30日收到修改稿)

福布斯下降(Forbush decrease,FD)是銀河宇宙線(galactic cosm ic rays,GCRs)受短期劇烈太陽活動(dòng)調(diào)制的重要現(xiàn)象之一.本文設(shè)GCRs進(jìn)入由行星際日冕物質(zhì)拋射(interplanetary coronalmass ejection,ICME)及其前沿激波共同形成的擾動(dòng)區(qū)時(shí),其徑向擴(kuò)散系數(shù)κrr受抑制變?yōu)棣?r)·κrr(0＜μ(r)≤1),且抑制強(qiáng)度與粒子位置處的太陽風(fēng)等離子體速度正相關(guān).對(duì)任意時(shí)刻的擾動(dòng)區(qū),抑制系數(shù)μ(r)在激波處最小為μ(rsh),并按指數(shù)規(guī)律增大,在ICME尾部歸一.CME爆發(fā)時(shí),μ(rsh)取全局最小值μm.在擾動(dòng)區(qū)向日球?qū)油鈧鞑サ倪^程中,μ(rsh)逐步恢復(fù)為1.在此基礎(chǔ)上,根據(jù)GOES和ACE衛(wèi)星觀測(cè)確定模型參數(shù),用一維隨機(jī)微分方程描述GCRs在日球?qū)觾?nèi)的傳播,并采用倒向隨機(jī)方法模擬了一個(gè)由獨(dú)立Halo ICME調(diào)制GCRs引起的2005年5月15日FD事件.計(jì)算所得地面中子通量的主相、恢復(fù)相及其在CME到達(dá)地球前的增加過程,均與Oulu中子探測(cè)器觀測(cè)結(jié)果一致.

行星際日冕物質(zhì)拋射,福布斯下降,倒向隨機(jī)微分方法,中子通量

1 引 言

銀河宇宙線(galactic cosm ic rays,GCRs)是起源于太陽系之外,主要由質(zhì)子、α粒子和少量電子組成的高能粒子,其能譜基本服從冪律分布,能量可達(dá)到1022eV[1,2].高能GCRs穿越由太陽風(fēng)等離子體形成的日球?qū)?到達(dá)地球大氣層后,與大氣發(fā)生碰撞并使氣體分子電離,形成廣延大氣簇射.研究顯示,GCRs是20 km以下大氣的主要電離源,能導(dǎo)致直接和間接的空間輻射事件,也與地磁強(qiáng)度等空間環(huán)境要素存在顯著關(guān)聯(lián)[3,4].

由太陽活動(dòng)爆發(fā)引起的不同時(shí)間尺度和強(qiáng)度的等離子體擾動(dòng),會(huì)影響GCRs在日球?qū)又械膫鞑?太陽活動(dòng)劇烈時(shí),到達(dá)地球的GCRs受到抑制,反之GCRs強(qiáng)度會(huì)增強(qiáng).太陽活動(dòng)對(duì)GCRs的調(diào)制效應(yīng)可根據(jù)調(diào)制因素和時(shí)間尺度分為很多種[5].其中福布斯下降(Forbush decrease,FD)事件是由短時(shí)太陽劇烈活動(dòng)導(dǎo)致地球位置GCRs通量急劇下降并逐漸恢復(fù)的現(xiàn)象.常見的FD事件為非重現(xiàn)型事件,下降相和恢復(fù)相不對(duì)稱,其通量下降相一般持續(xù)幾小時(shí)或一兩天,恢復(fù)相是隨后幾天或十幾天[6,7].

觀測(cè)顯示,FD事件與太陽日冕物質(zhì)拋射(coronal mass ejection,CME)引起的行星際激波和磁云關(guān)系密切[5,8-11].對(duì)FD事件中地球位置GCRs通量變化的研究,有助于人們?nèi)嬲J(rèn)識(shí)ICME對(duì)日地空間環(huán)境的影響,也可反推ICME相關(guān)的空間等離子體狀態(tài)變化[5,12].地面GCRs通量變化與地磁活動(dòng)Dst指數(shù)之間存在明顯關(guān)聯(lián),因此FD與磁暴事件都可以作為空間環(huán)境變化的預(yù)報(bào)指針[10,13].近地空間雷暴與FD事件也存在一定的相關(guān)性[14].弄清FD事件的根源,及其與太陽活動(dòng)、行星際磁場(chǎng)結(jié)構(gòu)、電離層狀態(tài)間的關(guān)系,有利于提高空間天氣和空間環(huán)境預(yù)測(cè)和預(yù)警水平.

為理解太陽活動(dòng)對(duì)GCRs的調(diào)制現(xiàn)象,人們提出了各類模型[5-8,15-24].早在1971年,Fisk[19]就建立了一維傳輸方程的穩(wěn)態(tài)數(shù)值解.隨著觀測(cè)水平的提高,包含多種行星際要素的計(jì)算模擬研究也越來越多[20-22].其中GCRs調(diào)制的力場(chǎng)模型用作用勢(shì)參數(shù)描述太陽風(fēng)對(duì)GCRs輸運(yùn)的調(diào)制,給出了GCRs粒子調(diào)制能譜變化的唯象描述[23-25].力場(chǎng)模型在GCRs短期調(diào)制中的應(yīng)用也已經(jīng)引起人們的興趣[26].

在GCRs日球?qū)虞斶\(yùn)模擬中,隨機(jī)微分方程(stochastic differential equation,SDE)方法因數(shù)值處理方便、物理圖像簡明[27,28],被越來越多的課題組采用[17,18,27-33].Waw rzynczak等[29]基于GCRs日球?qū)觽鬏敺匠?利用SDE方法建立了GCRs的FD現(xiàn)象和27天變化的短時(shí)調(diào)制模型,并將結(jié)果與有限差分方法比較,發(fā)現(xiàn)兩種方法都與觀測(cè)結(jié)果符合.Pei等[30]通過建立球坐標(biāo)系中的三維含時(shí)SDE方程,求解Parker傳播方程,提高了程序的有效性.Li等[17,18]曾嘗試通過一維隨機(jī)方程模型建立宇宙射線被太陽活動(dòng)調(diào)制的直觀物理機(jī)制. Luo等[31,32]基于SDE方法,研究了全局耦合作用對(duì)GCRs的調(diào)制效應(yīng).Bobik等[33]系統(tǒng)研究了第23太陽周中GCRs質(zhì)子的太陽調(diào)制,分析了漂移和緯度效應(yīng).

由于ICME等太陽活動(dòng)爆發(fā)過程的復(fù)雜性,細(xì)致研究FD事件與空間等離子體狀態(tài)變化的關(guān)聯(lián)仍比較困難.在空間天氣過程分析中,對(duì)FD事件產(chǎn)生和演化的理解,仍需要一個(gè)直觀、簡單可靠的調(diào)制模型.在CME與FD事件的基本關(guān)聯(lián)模型的基礎(chǔ)上[8],本文用一維SDE方法模擬GCRs在日球?qū)觾?nèi)的傳播,用GCRs擴(kuò)散系數(shù)κrr對(duì)ICME激波和磁云區(qū)的響應(yīng)函數(shù)描述單個(gè)ICME對(duì)GCRs的調(diào)制,計(jì)算了一個(gè)典型的非重現(xiàn)型FD事件發(fā)生過程中地面中子通量的變化,并將計(jì)算結(jié)果與中子探測(cè)器觀測(cè)進(jìn)行了比較和分析.

本文第2部分給出GCRs一維輸運(yùn)的SDE方程,并介紹ICME調(diào)制GCRs擴(kuò)散的基本模型;第3部分以GOES和ACE衛(wèi)星觀測(cè)為基礎(chǔ),通過設(shè)置合理的調(diào)制模型參數(shù)組,模擬并計(jì)算了2005年5月15日FD事件中地面中子通量的演化;最后,是本工作的總結(jié)與討論.

2 物理模型與模擬方法

2.1 GCRs傳播的隨機(jī)微分方程

Parker的宇宙射線傳輸方程為[34]

其中,f是GCRs的時(shí)空分布函數(shù),t是時(shí)間,Vsw是太陽風(fēng)速度,κ是GCRs的擴(kuò)散系數(shù),P是粒子的剛度(P=pc/q,p為粒子動(dòng)量大小,q為粒子所帶電荷).等號(hào)右邊第一項(xiàng)為對(duì)流項(xiàng),第二項(xiàng)為擴(kuò)散項(xiàng),第三項(xiàng)為絕熱能量改變項(xiàng).設(shè)方程(1)具有球?qū)ΨQ性,在以太陽中心為原點(diǎn)、日地連線為徑向的一維坐標(biāo)系中,方程(1)簡化為[17]

其中,r是粒子到太陽中心的距離,κrr是粒子的一維擴(kuò)散系數(shù).

通過求解上述擴(kuò)散方程,理論上可求得GCRs的宏觀分布函數(shù)f(t,r,p).這種宏觀平均雖然可直接觀測(cè),但在復(fù)雜初值和邊界條件下,其求解過程并不容易.GCRs在日球?qū)拥妮斶\(yùn),本質(zhì)上可看作高能帶電粒子在行星際磁場(chǎng)中的隨機(jī)行走.雖然單個(gè)粒子的隨機(jī)行為很難給出GCRs通量或能譜等觀測(cè)結(jié)果,但如果運(yùn)用馬爾可夫隨機(jī)過程理論對(duì)大量粒子進(jìn)行跟蹤和統(tǒng)計(jì),則不僅可以求得上述擴(kuò)散方程的解,還可以獲得諸如粒子軌跡、動(dòng)量損失、輸運(yùn)時(shí)間等更加詳細(xì)的信息,從而有效研究GCRs的調(diào)制效應(yīng)[28,32].設(shè)GCRs在日球?qū)油膺吔鐁o處具有確定的能譜,在ro處釋放試探粒子,跟蹤其在日球?qū)觾?nèi)的隨機(jī)行走至1 AU處,記錄任意時(shí)刻1 AU處粒子的數(shù)量和能量,即可統(tǒng)計(jì)出地球附近的GCRs通量.

為提高統(tǒng)計(jì)效率,實(shí)際可采用倒向隨機(jī)方法[28,31]：將1 AU處能量為E的粒子倒向回溯至日球?qū)油膺吔鐁o處,并記錄其邊界處能量E′.GCRs在日球?qū)油獾哪茏Vj(E′)正是能量為E′的粒子自日球?qū)油庖阅芰縀到達(dá)1 AU的概率.對(duì)1 AU處可能觀測(cè)到的不同能量粒子進(jìn)行回溯模擬,可求得

該點(diǎn)GCRs相對(duì)分布：

其中,N(E)是t時(shí)刻在1 AU處釋放的能量為E的粒子數(shù)總和.

在實(shí)際的倒向SDE中,令d s=-d t(d t＜0),則(2)式可改寫成描述準(zhǔn)粒子運(yùn)動(dòng)的隨機(jī)微分方程組[17,35]：

其中,ξ為均勻分布在(0,1)之間的隨機(jī)數(shù);erf-1(ξ)是余誤差函數(shù);表示維納過程,保證了粒子的隨機(jī)行走過程.上式推導(dǎo)時(shí),假設(shè)Vsw與徑向無關(guān)[17,36].忽略日鞘及終止激波的影響,可設(shè)日球?qū)油膺吔缣嶨CRs通量分布規(guī)律為：j(E′)～(E′)-2.35.

2.2 GCRs擴(kuò)散系數(shù)κrr對(duì)ICM E擾動(dòng)的響應(yīng)

設(shè)日球?qū)油釭CRs的擴(kuò)散行為具有各向同性.若無太陽活動(dòng),粒子徑向擴(kuò)散系數(shù)κrr不隨時(shí)間變化.設(shè)粒子的靜止能量為E0,速度為v,動(dòng)量為p,質(zhì)量為m,質(zhì)量數(shù)為A,光速為c,則其擴(kuò)散系數(shù)可表示為：κrr=6×1022β(P/1GV)cm2/s.其中,

在一維力場(chǎng)模型中,用等離子體作用勢(shì)參數(shù)描述太陽活動(dòng)對(duì)粒子動(dòng)能的抑制效應(yīng)[23,24].Li等[17,18]以及Chih和Lee[37]用GCRs擴(kuò)散系數(shù)κrr的余弦響應(yīng)函數(shù)表示11年周期性太陽調(diào)制,研究了O+8通量變化的滯后現(xiàn)象.本文通過設(shè)計(jì)擴(kuò)散系數(shù)κrr對(duì)太陽活動(dòng)擾動(dòng)的動(dòng)態(tài)響應(yīng)函數(shù),描述ICME對(duì)GCRs的調(diào)制作用,模擬FD事件的發(fā)生和發(fā)展過程.

如圖1(a)中ICME簡圖所示,CME爆發(fā)后,當(dāng)磁云速度大于背景太陽風(fēng)等離子體的快磁聲波速度時(shí),CME前方通常會(huì)形成一個(gè)由激波和磁鞘組成的湍流區(qū),CME磁云則可能通過磁力線與太陽表面相連[38,39].在日地連線上,由激波、磁鞘和磁云共同組成的GCRs擾動(dòng)區(qū)中,等離子體速度場(chǎng)結(jié)構(gòu)復(fù)雜,但其徑向分量V(r)的分布大致如圖1(b)所示：激波速度最高,其后磁云的速度逐漸下降,至CME尾部速度將接近環(huán)境太陽風(fēng)速度.高能GCRs粒子進(jìn)入激波和磁云,與等離子體波動(dòng)發(fā)生相互作用,可能發(fā)生復(fù)雜的粒子俘獲、加速、逃逸等物理過程[5,40].

圖1 (網(wǎng)刊彩色)ICME擾動(dòng)區(qū)結(jié)構(gòu)及其對(duì)GCRs擴(kuò)散抑制效應(yīng)示意圖 (a)ICME及其前向激波在日地空間傳播;(b)擾動(dòng)區(qū)各處太陽風(fēng)速度Vsw分布;(c)擾動(dòng)區(qū)對(duì)GCRs輸運(yùn)的抑制系數(shù)μ(r)Fig.1.(color on line)The schem es of the disturbed zone induced by ICM E and itsm odu lation to GCRs: (a)The schem e of ICME and the forward IP shock m oving in the solar-terrestrial space;(b)the radial solar wind speed Vsw of the GCRs barrier region;(c)the hold ing strengthμ(r)of GCRs transport.

設(shè)CME擾動(dòng)區(qū)對(duì)高能GCRs粒子的作用以輸運(yùn)抑制為主[5,17].擾動(dòng)區(qū)域?qū)ο鄬?duì)論粒子的俘獲效應(yīng)類似粒子蓄水池,使進(jìn)入該區(qū)的GCRs粒子擴(kuò)散能力明顯下降.用參數(shù)μ(r)(0＜μ(r)≤1)表示GCRs粒子擴(kuò)散水平在湍流區(qū)中的降低程度.粒子進(jìn)入激波湍流和磁云擾動(dòng)區(qū)后,其含時(shí)擴(kuò)散系數(shù)被調(diào)制為：為簡單起見,設(shè)抑制強(qiáng)度與粒子所處位置的太陽風(fēng)速度Vsw(r)正相關(guān),則μ(r)與Vsw(r)反相關(guān)(如圖1(c)所示).

對(duì)任意時(shí)刻的ICME擾動(dòng)區(qū)而言,激波和磁鞘區(qū)的太陽風(fēng)速度最大,對(duì)GCRs擴(kuò)散的抑制效應(yīng)最強(qiáng),因此抑制系數(shù)最小,用μ(rsh)表示,其中rsh是激波位置.隨著ICME向日球?qū)舆吔邕\(yùn)動(dòng),激波強(qiáng)度不斷衰減,μ(rsh)逐漸恢復(fù)至1.該過程可用開關(guān)函數(shù)表示為

其中,Rc為激波衰減截?cái)辔恢?當(dāng)激波越過Rc位置后,強(qiáng)度開始明顯減弱.m和n表示調(diào)制程度減弱的系數(shù),其數(shù)值由ICME爆發(fā)和傳播過程確定.在激波前沿(寬度為Wsh-front)區(qū)域,GCRs擴(kuò)散系數(shù)從平靜值κrr隨位置線性下降至μ(rsh)·κrr.在激波與磁鞘后方,從ICME前沿向磁云尾部方向, μ(r)從μ(rsh)逐漸增加,并在磁云尾部歸一,即GCRs擴(kuò)散系數(shù)恢復(fù)正常.假設(shè)ICME磁云中μ(r)的恢復(fù)過程為指數(shù)形式：μ(r)～exp(τ·rcme),其中rcme為粒子在CME磁云內(nèi)部的相對(duì)位置,τ為擴(kuò)散調(diào)制恢復(fù)指數(shù).

對(duì)于擾動(dòng)區(qū)寬度W,利用GOES和ACE衛(wèi)星觀測(cè),我們可確定1 AU處CME的寬度W1AU.又根據(jù)Wang等[41]的統(tǒng)計(jì)工作,在離太陽15 AU內(nèi),對(duì)地ICME的擾動(dòng)區(qū)寬度W隨其前沿激波位置rsh線性增加：W=0.05+0.16rsh.在15 AU外,可認(rèn)為擾動(dòng)區(qū)寬度基本保持不變.利用此CME寬度的平均變化率和1 AU處的觀測(cè)值,可確定ICME擾動(dòng)區(qū)從太陽出發(fā)逐步增寬并趨于恒定的過程.至此,我們建立了ICME調(diào)制GCRs的簡單參數(shù)化模型.該模型忽略粒子非徑向擴(kuò)散,只考慮擴(kuò)散系數(shù)的徑向變化;不考慮擾動(dòng)區(qū)三維結(jié)構(gòu),僅用簡單函數(shù)表示擾動(dòng)強(qiáng)度對(duì)ICME徑向結(jié)構(gòu)和傳播速度的響應(yīng).

3 FD事件模擬與討論

本文選取2005年5月15日FD事件為研究對(duì)象.觀測(cè)表明,該事件是由一個(gè)獨(dú)立的Halo CME引起[42,43].事件中ICME的角度分布較寬,對(duì)地有效性大,適合用一維模型來描述.另一方面,該FD事件獨(dú)立性較好,下降相和恢復(fù)相分區(qū)清晰,有利于對(duì)物理模型開展評(píng)估分析.

3.1 模型參數(shù)

圖2(1-5)分別給出了GOES衛(wèi)星觀測(cè)[44]到的X-ray強(qiáng)度和ACE衛(wèi)星觀測(cè)[45]到的太陽風(fēng)速度x分量Vx、磁場(chǎng)強(qiáng)度大小|B|及其z分量Bz,質(zhì)子數(shù)密度Np的演化過程.由X-ray的峰值位置可以看出,CME在5月13日15：12爆發(fā),伴隨耀斑為M級(jí). 5月15日2：18激波到達(dá)1 AU.CME在當(dāng)天6：00到 22：00間穿越1 AU位置.激波和CME先后到達(dá)1 AU的時(shí)差為3.7 h,此即激波磁鞘區(qū)域通過1 AU的時(shí)長.觀測(cè)表明,CME爆發(fā)速度為1240.5 km/s,到達(dá)1 AU處速度為857 km/s,背景太陽風(fēng)速度為500 km/s.圖2(6)給出了Oulu中子探測(cè)器[46]記錄的地面中子通量相對(duì)變化.CME爆發(fā)后自激波到達(dá)1 AU起,地面中子通量在6 h內(nèi)下降幅度達(dá)到12%,并在隨后十天逐漸恢復(fù).因此,該Halo ICME引發(fā)了一個(gè)典型的非重現(xiàn)型FD事件.

表1 GCRs擴(kuò)散調(diào)制模型中的參數(shù)物理意義及其在2005年5月15日FD事件中的數(shù)值或相對(duì)大小Tab le 1.The param eters of the GCRs diff usion m odu lation model,and their values for the FD event on May 15,2005.

由于GCRs通量主要由高能質(zhì)子確定,本文以質(zhì)子為試探粒子,模擬ICME引起FD事件的發(fā)生過程.模擬中設(shè)日球?qū)觾?nèi)邊界rin=0.02 AU,外邊界ro=95 AU.模擬t時(shí)刻通量時(shí),在觀測(cè)點(diǎn)(1 AU)處釋放大量的準(zhǔn)粒子,使其倒向隨機(jī)擴(kuò)散直到外邊界ro外,運(yùn)動(dòng)軌跡遵循隨機(jī)微分方程組(3)和(4).若粒子運(yùn)動(dòng)到內(nèi)邊界rin內(nèi),對(duì)粒子進(jìn)行鏡像對(duì)稱操作,使其返回到日球?qū)觾?nèi)繼續(xù)運(yùn)動(dòng).

圖2 (網(wǎng)刊彩色)2005年5月13日到2005年5月16日期間GOES衛(wèi)星觀測(cè)到的X-ray強(qiáng)度和ACE衛(wèi)星觀測(cè)到的太陽風(fēng)速度分量Vx、磁場(chǎng)強(qiáng)度|B|及其z分量Bz,質(zhì)子數(shù)密度Np,以及Ou lu中子探測(cè)器記錄的地面中子通量的演化Fig.2.(color on line)The X-ray intensity recorded by GOES, the rad ial speed of solar wind Vx,the m agnetic fi led|B| and its z com ponent Bz,the proton density Np recorded by ACE,and the neu tron fl ux at the ground level recorded by the Ou lu neu tron m onitor,from May 13 to 15 of 2005.

模擬中準(zhǔn)粒子能量取值范圍以O(shè)ulu中子探測(cè)器臺(tái)站截止剛度為參考,取為300 MeV-150 GeV.按照能量對(duì)數(shù)等間距原則,分成20個(gè)測(cè)試點(diǎn).每個(gè)能量點(diǎn)在同一觀測(cè)時(shí)刻釋放足夠多的準(zhǔn)粒子,確保隨機(jī)結(jié)果漲落與觀測(cè)基本一致.地面中子通量N由1 AU處初級(jí)宇宙線譜與產(chǎn)額函數(shù)的卷積決定[47,48]：

其中,Pc表示當(dāng)?shù)氐卮沤刂箘偠?h是大氣深度; Ji(P,t)[GV m2·sr·s]是初級(jí)宇宙線i在時(shí)刻t的剛度譜;Yi(P,h)[m2·sr]是初級(jí)宇宙線i對(duì)應(yīng)的中子探測(cè)器的產(chǎn)額函數(shù).其表達(dá)式為[47,48]：

其中,Ai(E,θ)是探測(cè)面積與計(jì)數(shù)率之積;Fi,j表示次級(jí)宇宙線粒子j(有中子、質(zhì)子、介子等)的微分通量,θ是次級(jí)宇宙線的入射角.本文采用文獻(xiàn)[47,48]報(bào)道的Oulu臺(tái)站的質(zhì)子-中子參數(shù)化產(chǎn)額函數(shù),由模擬所得1 AU處GCRs通量,計(jì)算出中子通量的相對(duì)變化,并與觀測(cè)進(jìn)行對(duì)比.

3.2 模擬與計(jì)算結(jié)果

圖3(a)中藍(lán)色實(shí)線給出了幾個(gè)等間距觀測(cè)時(shí)刻,1 GeV質(zhì)子在1-5 AU內(nèi)的倒向行走軌跡以及激波磁鞘(紅色),CME前沿和磁云尾部半寬(綠色)的時(shí)間演化過程,可看出GCRs的倒向隨機(jī)擴(kuò)散行為和擾動(dòng)區(qū)形態(tài)演化特征.圖3(b)給出了根據(jù)模擬所得1 AU處的GCRs通量,經(jīng)Oulu臺(tái)站中子產(chǎn)額函數(shù)變換[47,48],得到的地面中子通量(紅色)與Oulu中子探測(cè)器的15m in分辨率記錄值(黑色).GCRs粒子通量在5月15日01時(shí)開始下降,地面中子通量亦隨之降低,在5月15日09時(shí)降到最低.其后,中子通量開始緩慢恢復(fù).恢復(fù)相從5月15日09時(shí)到5月20日12時(shí)大約跨越5天零3小時(shí).計(jì)算結(jié)果與實(shí)際觀測(cè)相符.

從圖3(a)中試探粒子的軌跡看出,在ICME傳播至1 AU之前,到達(dá)地球的GCRs粒子幾乎未受ICME影響,因此中子通量保持穩(wěn)定.由于前向激波的前沿具有擋板效應(yīng),在ICME到達(dá)1 AU處前,中子通量急劇下降前出現(xiàn)了約1%的前期增加[49].測(cè)試表明,本文模型計(jì)算所得中子通量的前期增幅與激波前沿寬度Wsh-front正相關(guān).在ICME通過1 AU后,到達(dá)地球的GCRs粒子均經(jīng)歷了ICME擾動(dòng)區(qū)域的調(diào)制.當(dāng)擾動(dòng)區(qū)運(yùn)動(dòng)到1-3 AU區(qū)域時(shí),磁鞘區(qū)中粒子滯留時(shí)間相對(duì)較長,粒子隨磁鞘湍流漂移導(dǎo)致粒子到達(dá)地球的時(shí)間被顯著推移.隨著ICME的傳播,激波強(qiáng)度逐漸減弱,且CME寬度逐漸增加,粒子在磁云中滯留(往復(fù)運(yùn)動(dòng))的時(shí)間有所增加,在圖3(a)中表現(xiàn)為粒子在磁云區(qū)域的振蕩次數(shù)的增多.隨著擾動(dòng)區(qū)GCRs粒子數(shù)目變多,而擾動(dòng)強(qiáng)度卻不斷衰減,出入ICME區(qū)域的粒子數(shù)趨向平衡,地面中子通量亦隨之恢復(fù)至寧靜水平.

圖3 (網(wǎng)刊彩色)(a)2005年5月15日FD事件模擬中1 GeV粒子的輸運(yùn)軌跡與ICME調(diào)制區(qū)傳播過程;(b)計(jì)算所得中子通量與Ou lu臺(tái)站觀測(cè)對(duì)比Fig.3.(color on line)(a)The tra jectroise of 1 GeV test particles and the evolution of the ICME disturbing area du ring ou r sim u lation of FD event on the 15 May,2005;(b)the com parison d iagram of the neu tron fl ux between Ou lu observation and ou r sim u lation.

從上述過程可進(jìn)一步理解本模型中的相關(guān)參數(shù)的取值.在FD事件發(fā)生期間,ICME擾動(dòng)區(qū)在行星際移動(dòng)的有效距離約3 AU,這與模型中ICME調(diào)制強(qiáng)度的截?cái)嗑嚯xRc相近.當(dāng)擾動(dòng)區(qū)運(yùn)動(dòng)到3 AU之外后,ICME對(duì)GCRs的俘獲能力逐漸達(dá)到飽和并開始衰減,因此用開關(guān)函數(shù)表達(dá)激波對(duì)GCRs擴(kuò)散系數(shù)的調(diào)制效應(yīng)具有一定的合理性.地面中子通量曲線的快速下降和緩慢恢復(fù)的趨勢(shì),實(shí)際上反映了太陽活動(dòng)的調(diào)制區(qū)域?qū)CRs輸運(yùn)的抑制效應(yīng)先強(qiáng)后弱的動(dòng)態(tài)特性.

4 總結(jié)與討論

FD事件是太陽活動(dòng)爆發(fā)調(diào)制GCRs傳播,導(dǎo)致1 AU處GCRs和其地面次級(jí)中子通量變化的空間天氣過程.本文用GCRs在行星際日地徑向擴(kuò)散系數(shù)κrr的變化,表示ICME引發(fā)的空間等離子體擾動(dòng)對(duì)GCRs的調(diào)制.高能粒子的徑向擴(kuò)散抑制,本質(zhì)上可導(dǎo)致其沿太陽風(fēng)向日球?qū)油舛ㄏ蚱?設(shè)調(diào)制強(qiáng)度與擾動(dòng)區(qū)太陽風(fēng)徑向速度正相關(guān).在擾動(dòng)區(qū)向日球?qū)油鈧鞑ミ^程中,調(diào)制效應(yīng)不斷減弱.將ICME前向激波對(duì)GCRs的調(diào)制強(qiáng)度用開關(guān)函數(shù)表示;在磁云區(qū)內(nèi)部,自CME前沿至其尾部,調(diào)制強(qiáng)度呈指數(shù)衰減.在此參數(shù)化模型的基礎(chǔ)上,用一維隨機(jī)微分方程模擬1 AU處的GCRs通量變化,并計(jì)算地面中子通量曲線.對(duì)2005年5月13日爆發(fā)的單Halo CME引起的5月15日FD事件,該模型計(jì)算得到的地面中子通量與Oulu中子探測(cè)器記錄的變化趨勢(shì)一致,成功模擬了FD事件發(fā)生前中子通量的前期增加、激波經(jīng)過1 AU處時(shí)中子通量的迅速下降(主相)和其后的逐漸恢復(fù)過程(恢復(fù)相).

在本文ICME導(dǎo)致FD事件的一維模型中,太陽風(fēng)和ICME結(jié)構(gòu)的確定主要以GOES和ACE衛(wèi)星測(cè)得的CME爆發(fā)及其通過1 AU處時(shí)的等離子體和磁場(chǎng)參數(shù)為依據(jù),調(diào)制強(qiáng)度μ(r)的變化需要適當(dāng)調(diào)節(jié).實(shí)際上,很多引起FD事件的CME未必具有較好的對(duì)地性,ICME結(jié)構(gòu)參數(shù)比較復(fù)雜[11].隨著多衛(wèi)星觀測(cè)和CME三維重構(gòu)技術(shù)的提升[50],可通過觀測(cè)和重構(gòu)獲得更加全面的ICME擾動(dòng)區(qū)調(diào)制信息.另一方面,模擬地面中子通量變化時(shí)對(duì)擾動(dòng)區(qū)參數(shù)的調(diào)節(jié),亦可看作是根據(jù)地面中子監(jiān)測(cè)對(duì)ICME結(jié)構(gòu)和強(qiáng)度變化進(jìn)行的反演和重構(gòu).以反演所得的ICME參數(shù)為依據(jù),我們可以模擬出其他位置(如R＞1 AU,火星位置)的GCRs通量變化.這對(duì)行星際空間天氣預(yù)報(bào)有一定的參考價(jià)值.

為進(jìn)一步對(duì)FD通量曲線中的細(xì)節(jié)振蕩特征進(jìn)行研究,需要對(duì)太陽活動(dòng)調(diào)制GCRs過程進(jìn)行更加詳細(xì)的刻畫、采用包含更多物理因素的多維模型、考慮日球?qū)幼陨砗托行请H等離子體的三維結(jié)構(gòu).將SDE求解GCRs的傳播過程擴(kuò)展到三維空間,加入太陽風(fēng)磁場(chǎng)結(jié)構(gòu),以及粒子在ICME激波和磁云中的具體相互作用,研究行星際激波湍動(dòng)、電流片等因素對(duì)GCRs傳播的影響,也是我們努力的目標(biāo).

感謝美國阿拉巴馬大學(xué)(亨茨維爾)李剛教授給予的指導(dǎo)和幫助,感謝南京信息工程大學(xué)丁留貫副教授和紫金山天文臺(tái)封莉研究員的討論和幫助.

[1]Rossi B 1964 Cosm ic Rays(New York:M cG raw-H ill) pp1-10

[2]B lasi P 2013 Astron.Astrophys.Rev.21 70

[3]Bothm er V,Daglis IA 2007 Space W eather-Physics and Effects(Berlin:Springer)pp103-130

[4]Le G M 2002 Ph.D.D issertation(Beijing:Chinese Academ y of Sciences)(in Chinese)[樂貴明 2002博士學(xué)位論文(北京:中國科學(xué)院空間科學(xué)與應(yīng)用研究中心)]

[5]Potgieter M S 2013 Living Rev.Sol.Phys.10 3

[6]Guo W J,Zhu B Y 1990 Chin.J.Spac.Sci.10 247(in Chinese)[郭維吉,朱邦耀1990空間科學(xué)學(xué)報(bào)10 247]

[7]K harayat H,Prasad L,M athpal R,Garia S,Bhatt B 2016 Sol.Phys.291 603

[8]Cane H V 2000 Space Sci.Rev.93 55

[9]Belov A,Abunin A,Abunina M,Eroshenko E,O leneva V,Yanke V,Papaioannou A,M avrom ichalaki H,Gopalswam y N,Yashiro S 2014 So l.Phys.289 3949

[10]Yu X X,Lu H,Le G M,Shi F 2010 So l.Phys.263 223

[11]Zhao L L,Zhang H 2016 Astrophys.J.827 13

[12]Le G M,Han Y B 2005 Acta Phys.Sin.54 467(in Chinese)[樂貴明,韓延本2005物理學(xué)報(bào)54 467]

[13]Dorm an L I 2005 Ann.Geophys.23 2997

[14]Huang Y L,Fu Y F,Chen J M,Huang G S,Liu X N 2015 Res.Astron.Astrophys.16 82(in Chinese)[黃寅亮,傅元芬,陳濟(jì)民,黃更生,劉小寧2015天體物理學(xué)報(bào)16 82]

[15]Lockwood J A 1971 Space Sci.Rev.12 658

[16]Jokip ii J R,Kota J 1986 J.Geophys.Res.91 2885

[17]Li G,W ebb G M,Roux J A L,Zank G P,W iedenbeck M E 2007 Num erical M odeling of Space Plasm a Flows 385 31

[18]LiG,W ebb G M,Roux J A L,W iedenbeck M,F lorinski V,Zank G P 2009 Proceedings of the 31st ICRC Tód?, Poland,Ju ly 7-15,2009

[19]Fisk L A 1971 J.Geophys.Res.76 221

[20]Hattingh M,Burger R A,Potgieter M S,Haasb roeket L J 1997 Adv.Space Res.19 893

[21]G il A,Iskra K,M odzelew ska R,Alania M V 2005 Adv. Space Res.35 687

[22]Strauss R D,Potgieter M S,Büsching I,Kopp A 2012 Astrophys.Space Sci.339 223

[23]G leeson L 1968 Astrophys.J.154 1011

[24]Caballero-Lopez R A,M oraal H,M cdonald F B 2004 J. Geophys.Res.109 361

[25]Usoskin IG,Bazilevskaya G A,Kovaltsov G A 2011 J. Geophys.Res.116 1

[26]Usoskin IG,Kovaltsov G A,AdrianiO,Barbarino G C, Bazilevskaya G A,Bellotti R 2015 Adv.Space Res.55 2940

[27]Dunzlaff P,Strauss R D,Potgieter M S 2015 Com put. Phys.Comm un.192 156

[28]Zhang M 1999 Astrophys.J.513 409

[29]W aw rzynczak A,M odzelew ska R,Gil A 2015 J.Phys. Conf.Ser.574 012078

[30]Pei C,Bieber J W,B reech B,Bu rger R A,C lem J, M atthaeus W H 2010 J.Geophys.Res.115 333

[31]Luo X,Zhang M,Rassou l H K,Pogorelov N V 2011 Astrophys.J.730 13

[32]Luo X,Zhang M,Feng X,M endoza-Torres J E 2013 J. Geophys.Res.118 7517

[33]Bobik P,Boella G,Boschini M J,Consoland i C,Della Torre S,Gervasi M,G rand i D,Kudela K,Pensotti S, Rancoita P G,TacconiM 2012 Astrophys.J.745 132

[34]Parker E N 1965 Planet.Space Sci.13 9

[35]Schuss Z 1980 Theory and Applications of Stochastic D iff eren tial Equations(New York:John W iley)pp10-95

[36]W iedenbeck M E,Davis A J,Leske R A,Binna W R, Cohen C M S,Cumm ings A C,de Nolfo G,Israel M H, Lab rador A W,M ewald t R A,Scott L M,Stone E C, von Rosenvinge T T 2005 Proceedings of the 29th International Cosm ic Ray Conference Pune,Ind ia,August 3-10,2005 p227

[37]Chih P P,Lee M A 1986 J.Geophys.Res.91 2903

[38]N ishida A 1982 J.Geophys.Res.87 6003

[39]Richardson IG,Cane H V 2010 Sol.Phys.264 189

[40]Thom as B T,Gall R 1984 J.Geophys.Res.89 2991

[41]W ang C,Du D,Richardson J D,Liu Y 2005 Proceedings of the Solar W ind 11/SOHO16,“Connecting Sun and Heliosphere”Conference(ESA SP-592)W histler, Canada,June 12-17,2005 p781

[42]Verm a P L,Patel N K,Pra japatiM 2014 J.Phys.Conf. Ser.511 012057

[43]Ah luwalia H S,Alania M V,W aw rzynczak A,Ygbuhay R C,FikaniM M 2014 Sol.Phys.289 1763

[44]http://sp id rngdcnoaagov/sp id r/datasetdo[2017-3-6]

[45]http://cdawgsfcnasagov/CM E l ist/[2017-3-6]

[46]http://cosm icraysou lu fi/[2017-3-6]

[47]M ishev A L,Usoskin IG,Kovalstov G A 2013 J.Geophys.Res.118 2783

[48]M angeard P S,Ru ff olo D,Saiz A,M ad lee S,Nu taro T 2016 J.Geophys.Res.121 7435

[49]Papailiou M,M av rom ichalaki H,Belov A,Eroshenko E, Yanke V 2012 Sol.Phys.276 337

[50]Feng L,Inhester B,W ei Y,Gan W Q,Zhang T L,W ang M Y 2012 Astrophys.J.751 18

Interplanetary coronal mass ejection induced forbush decrease event: a simulation study with one-dimensional stochastic differential method

Ni Su-Lan1)2)Gu Bin1)2)?Han Zhi-Yi1)2)

1)(Departm ent of Physics,Nanjing University of Inform ation Science and Technology,Nanjing 210044,China)
2)(Institute of Space W eather,Nanjing University of Inform ation Science and Technology,Nanjing 210044,China)
(Received 6 March 2017;revised manuscript received 30 March 2017)

Forbush decrease(FD)event is one of themost im portant short-term modulations of galactic cosm ic rays(GCRs) caused by intense solar activities such as interp lanetary coronalm ass ejection(ICM E).Them odulation m echanism s of GCRs by the disturbed interp lanetary magnetic fields(IM F)of ICME and the accom panying forward interp lanetary shock(IP)are not clear yet.

In this work,we present a one-dim ensional dynam icmodel of the GCR barrier driven by ICME.In our model,the time dependent radial diff usion coefficientκrrof GCRs is dep ressed to beμ(r)·κrr(0＜μ(r)≤1)as they run into the disturbed IMF.The scale factorμ(r)is inversely proportional to the local solar wind speed away from the Sun.W ithin the disturbed area at any time,μ(r)increases exponentially from the localm inimumμ(rsh)at the IP front to 1 at the end of the ICME tail.In addition,μ(rsh)sw itches gradually from its globalm inimumμmat the bursting of the CME to 1 as the shock m oving toward the outer boundary of the heliosphere.The geom etrical and dynam ic param eters of the ICME and IP are derived from the observations of GOES and ACE satellites.

Based on the stochastic transport theory,the one-dimensional backward stochastic differential equation(SDE) method is adopted to simulate the transport of GCRsm odulated by single halo ICME.The evolution of the neutron flux at the ground is calcu lated according to the recently reported proton-neutron yield function.As an exam p le,the FD event on 15 May 2005,caused by the CME event bursting on 13 May 2005,is studied and simulated.The results show that the calcu lated neutron flux evolution,including not only the m ain and recovery phases,but also the preenhancement before the arriving of the CME at the Earth,is consistent with the observation of Oulu neutron monitor.

According to the trajectories of GCRs,it can be found that,the per-enhancem ent of the neutron fl ux is a resu lt of the scattering by the forward IP passing 1 AU.Before the IP reaches the sw itch cutoff Rc,GCRs are evidently con fined in the sheath between the IP and CME.A fter that,the GCRs w ill stay for longer time in themagnetic cloud of the ICME as a result of the dam ping of IP strength.

The param eterzed one-dim ensional GCRsm odulation model and the SDE method,as have been confi rm ed by the neutron monitor observation on the Earth,can be used further to calcu late and predict the GCRs fluxes of other p laces, such as the M ars,in the heliosphere.

interp lanetary coronal m ass ejection,Forbush decrease,backward stochastic differential method,neutron flux

PACS:96.50.S-,96.50.sh,96.60.ph,02.50.Ey DO I:10.7498/aps.66.139601

?通信作者.E-m ail:gubin@nuist.edu.cn

PACS:96.50.S-,96.50.sh,96.60.ph,02.50.Ey DO I:10.7498/aps.66.139601

?Corresponding author.E-m ail:gubin@nuist.edu.cn