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

倪素蘭顧斌韓智伊

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

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

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

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

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

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

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

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

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

1 引 言

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

由太陽活動爆發引起的不同時間尺度和強度的等離子體擾動,會影響GCRs在日球層中的傳播.太陽活動劇烈時,到達地球的GCRs受到抑制,反之GCRs強度會增強.太陽活動對GCRs的調制效應可根據調制因素和時間尺度分為很多種[5].其中福布斯下降(Forbush decrease,FD)事件是由短時太陽劇烈活動導致地球位置GCRs通量急劇下降并逐漸恢復的現象.常見的FD事件為非重現型事件,下降相和恢復相不對稱,其通量下降相一般持續幾小時或一兩天,恢復相是隨后幾天或十幾天[6,7].

觀測顯示,FD事件與太陽日冕物質拋射(coronal mass ejection,CME)引起的行星際激波和磁云關系密切[5,8-11].對FD事件中地球位置GCRs通量變化的研究,有助于人們全面認識ICME對日地空間環境的影響,也可反推ICME相關的空間等離子體狀態變化[5,12].地面GCRs通量變化與地磁活動Dst指數之間存在明顯關聯,因此FD與磁暴事件都可以作為空間環境變化的預報指針[10,13].近地空間雷暴與FD事件也存在一定的相關性[14].弄清FD事件的根源,及其與太陽活動、行星際磁場結構、電離層狀態間的關系,有利于提高空間天氣和空間環境預測和預警水平.

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

在GCRs日球層輸運模擬中,隨機微分方程(stochastic differential equation,SDE)方法因數值處理方便、物理圖像簡明[27,28],被越來越多的課題組采用[17,18,27-33].Waw rzynczak等[29]基于GCRs日球層傳輸方程,利用SDE方法建立了GCRs的FD現象和27天變化的短時調制模型,并將結果與有限差分方法比較,發現兩種方法都與觀測結果符合.Pei等[30]通過建立球坐標系中的三維含時SDE方程,求解Parker傳播方程,提高了程序的有效性.Li等[17,18]曾嘗試通過一維隨機方程模型建立宇宙射線被太陽活動調制的直觀物理機制. Luo等[31,32]基于SDE方法,研究了全局耦合作用對GCRs的調制效應.Bobik等[33]系統研究了第23太陽周中GCRs質子的太陽調制,分析了漂移和緯度效應.

由于ICME等太陽活動爆發過程的復雜性,細致研究FD事件與空間等離子體狀態變化的關聯仍比較困難.在空間天氣過程分析中,對FD事件產生和演化的理解,仍需要一個直觀、簡單可靠的調制模型.在CME與FD事件的基本關聯模型的基礎上[8],本文用一維SDE方法模擬GCRs在日球層內的傳播,用GCRs擴散系數κrr對ICME激波和磁云區的響應函數描述單個ICME對GCRs的調制,計算了一個典型的非重現型FD事件發生過程中地面中子通量的變化,并將計算結果與中子探測器觀測進行了比較和分析.

本文第2部分給出GCRs一維輸運的SDE方程,并介紹ICME調制GCRs擴散的基本模型;第3部分以GOES和ACE衛星觀測為基礎,通過設置合理的調制模型參數組,模擬并計算了2005年5月15日FD事件中地面中子通量的演化;最后,是本工作的總結與討論.

2 物理模型與模擬方法

2.1 GCRs傳播的隨機微分方程

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

其中,f是GCRs的時空分布函數,t是時間,Vsw是太陽風速度,κ是GCRs的擴散系數,P是粒子的剛度(P=pc/q,p為粒子動量大小,q為粒子所帶電荷).等號右邊第一項為對流項,第二項為擴散項,第三項為絕熱能量改變項.設方程(1)具有球對稱性,在以太陽中心為原點、日地連線為徑向的一維坐標系中,方程(1)簡化為[17]

其中,r是粒子到太陽中心的距離,κrr是粒子的一維擴散系數.

通過求解上述擴散方程,理論上可求得GCRs的宏觀分布函數f(t,r,p).這種宏觀平均雖然可直接觀測,但在復雜初值和邊界條件下,其求解過程并不容易.GCRs在日球層的輸運,本質上可看作高能帶電粒子在行星際磁場中的隨機行走.雖然單個粒子的隨機行為很難給出GCRs通量或能譜等觀測結果,但如果運用馬爾可夫隨機過程理論對大量粒子進行跟蹤和統計,則不僅可以求得上述擴散方程的解,還可以獲得諸如粒子軌跡、動量損失、輸運時間等更加詳細的信息,從而有效研究GCRs的調制效應[28,32].設GCRs在日球層外邊界ro處具有確定的能譜,在ro處釋放試探粒子,跟蹤其在日球層內的隨機行走至1 AU處,記錄任意時刻1 AU處粒子的數量和能量,即可統計出地球附近的GCRs通量.

為提高統計效率,實際可采用倒向隨機方法[28,31]：將1 AU處能量為E的粒子倒向回溯至日球層外邊界ro處,并記錄其邊界處能量E′.GCRs在日球層外的能譜j(E′)正是能量為E′的粒子自日球層外以能量E到達1 AU的概率.對1 AU處可能觀測到的不同能量粒子進行回溯模擬,可求得

該點GCRs相對分布：

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

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

其中,ξ為均勻分布在(0,1)之間的隨機數;erf-1(ξ)是余誤差函數;表示維納過程,保證了粒子的隨機行走過程.上式推導時,假設Vsw與徑向無關[17,36].忽略日鞘及終止激波的影響,可設日球層外邊界處GCRs通量分布規律為：j(E′)～(E′)-2.35.

2.2 GCRs擴散系數κrr對ICM E擾動的響應

設日球層外GCRs的擴散行為具有各向同性.若無太陽活動,粒子徑向擴散系數κrr不隨時間變化.設粒子的靜止能量為E0,速度為v,動量為p,質量為m,質量數為A,光速為c,則其擴散系數可表示為：κrr=6×1022β(P/1GV)cm2/s.其中,

在一維力場模型中,用等離子體作用勢參數描述太陽活動對粒子動能的抑制效應[23,24].Li等[17,18]以及Chih和Lee[37]用GCRs擴散系數κrr的余弦響應函數表示11年周期性太陽調制,研究了O+8通量變化的滯后現象.本文通過設計擴散系數κrr對太陽活動擾動的動態響應函數,描述ICME對GCRs的調制作用,模擬FD事件的發生和發展過程.

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

圖1 (網刊彩色)ICME擾動區結構及其對GCRs擴散抑制效應示意圖 (a)ICME及其前向激波在日地空間傳播;(b)擾動區各處太陽風速度Vsw分布;(c)擾動區對GCRs輸運的抑制系數μ(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.

設CME擾動區對高能GCRs粒子的作用以輸運抑制為主[5,17].擾動區域對相對論粒子的俘獲效應類似粒子蓄水池,使進入該區的GCRs粒子擴散能力明顯下降.用參數μ(r)(0＜μ(r)≤1)表示GCRs粒子擴散水平在湍流區中的降低程度.粒子進入激波湍流和磁云擾動區后,其含時擴散系數被調制為：為簡單起見,設抑制強度與粒子所處位置的太陽風速度Vsw(r)正相關,則μ(r)與Vsw(r)反相關(如圖1(c)所示).

對任意時刻的ICME擾動區而言,激波和磁鞘區的太陽風速度最大,對GCRs擴散的抑制效應最強,因此抑制系數最小,用μ(rsh)表示,其中rsh是激波位置.隨著ICME向日球層邊界運動,激波強度不斷衰減,μ(rsh)逐漸恢復至1.該過程可用開關函數表示為

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

對于擾動區寬度W,利用GOES和ACE衛星觀測,我們可確定1 AU處CME的寬度W1AU.又根據Wang等[41]的統計工作,在離太陽15 AU內,對地ICME的擾動區寬度W隨其前沿激波位置rsh線性增加：W=0.05+0.16rsh.在15 AU外,可認為擾動區寬度基本保持不變.利用此CME寬度的平均變化率和1 AU處的觀測值,可確定ICME擾動區從太陽出發逐步增寬并趨于恒定的過程.至此,我們建立了ICME調制GCRs的簡單參數化模型.該模型忽略粒子非徑向擴散,只考慮擴散系數的徑向變化;不考慮擾動區三維結構,僅用簡單函數表示擾動強度對ICME徑向結構和傳播速度的響應.

3 FD事件模擬與討論

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

3.1 模型參數

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

表1 GCRs擴散調制模型中的參數物理意義及其在2005年5月15日FD事件中的數值或相對大小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通量主要由高能質子確定,本文以質子為試探粒子,模擬ICME引起FD事件的發生過程.模擬中設日球層內邊界rin=0.02 AU,外邊界ro=95 AU.模擬t時刻通量時,在觀測點(1 AU)處釋放大量的準粒子,使其倒向隨機擴散直到外邊界ro外,運動軌跡遵循隨機微分方程組(3)和(4).若粒子運動到內邊界rin內,對粒子進行鏡像對稱操作,使其返回到日球層內繼續運動.

圖2 (網刊彩色)2005年5月13日到2005年5月16日期間GOES衛星觀測到的X-ray強度和ACE衛星觀測到的太陽風速度分量Vx、磁場強度|B|及其z分量Bz,質子數密度Np,以及Ou lu中子探測器記錄的地面中子通量的演化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.

模擬中準粒子能量取值范圍以Oulu中子探測器臺站截止剛度為參考,取為300 MeV-150 GeV.按照能量對數等間距原則,分成20個測試點.每個能量點在同一觀測時刻釋放足夠多的準粒子,確保隨機結果漲落與觀測基本一致.地面中子通量N由1 AU處初級宇宙線譜與產額函數的卷積決定[47,48]：

其中,Pc表示當地地磁截止剛度;h是大氣深度; Ji(P,t)[GV m2·sr·s]是初級宇宙線i在時刻t的剛度譜;Yi(P,h)[m2·sr]是初級宇宙線i對應的中子探測器的產額函數.其表達式為[47,48]：

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

3.2 模擬與計算結果

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

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

圖3 (網刊彩色)(a)2005年5月15日FD事件模擬中1 GeV粒子的輸運軌跡與ICME調制區傳播過程;(b)計算所得中子通量與Ou lu臺站觀測對比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.

從上述過程可進一步理解本模型中的相關參數的取值.在FD事件發生期間,ICME擾動區在行星際移動的有效距離約3 AU,這與模型中ICME調制強度的截斷距離Rc相近.當擾動區運動到3 AU之外后,ICME對GCRs的俘獲能力逐漸達到飽和并開始衰減,因此用開關函數表達激波對GCRs擴散系數的調制效應具有一定的合理性.地面中子通量曲線的快速下降和緩慢恢復的趨勢,實際上反映了太陽活動的調制區域對GCRs輸運的抑制效應先強后弱的動態特性.

4 總結與討論

FD事件是太陽活動爆發調制GCRs傳播,導致1 AU處GCRs和其地面次級中子通量變化的空間天氣過程.本文用GCRs在行星際日地徑向擴散系數κrr的變化,表示ICME引發的空間等離子體擾動對GCRs的調制.高能粒子的徑向擴散抑制,本質上可導致其沿太陽風向日球層外定向漂移.設調制強度與擾動區太陽風徑向速度正相關.在擾動區向日球層外傳播過程中,調制效應不斷減弱.將ICME前向激波對GCRs的調制強度用開關函數表示;在磁云區內部,自CME前沿至其尾部,調制強度呈指數衰減.在此參數化模型的基礎上,用一維隨機微分方程模擬1 AU處的GCRs通量變化,并計算地面中子通量曲線.對2005年5月13日爆發的單Halo CME引起的5月15日FD事件,該模型計算得到的地面中子通量與Oulu中子探測器記錄的變化趨勢一致,成功模擬了FD事件發生前中子通量的前期增加、激波經過1 AU處時中子通量的迅速下降(主相)和其后的逐漸恢復過程(恢復相).

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

為進一步對FD通量曲線中的細節振蕩特征進行研究,需要對太陽活動調制GCRs過程進行更加詳細的刻畫、采用包含更多物理因素的多維模型、考慮日球層自身和行星際等離子體的三維結構.將SDE求解GCRs的傳播過程擴展到三維空間,加入太陽風磁場結構,以及粒子在ICME激波和磁云中的具體相互作用,研究行星際激波湍動、電流片等因素對GCRs傳播的影響,也是我們努力的目標.

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

[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博士學位論文(北京:中國科學院空間科學與應用研究中心)]

[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空間科學學報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物理學報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)[黃寅亮,傅元芬,陳濟民,黃更生,劉小寧2015天體物理學報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