行星际日冕物质抛射引起福布斯下降的一维随机微分模拟

倪素兰顾斌韩智伊

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