内辐射带高能质子对磁暴响应的统计研究

何兆海，徐寄遥,2，王赤，戴磊，陈涛，Ilan Roth

1 中国科学院国家空间科学中心空间天气学国家重点实验室，北京 100080 2 中国科学院大学，北京 100049 3 Space Sciences Laboratory, University of California, Berkeley, California, USA

0 引言

早先研究普遍认为地球内辐射带相当稳定.然而，近来的观测证实磁暴和太阳高能粒子事件期间，内辐射带外边界的高能质子通量变化显著(Engel et al.，2015，2016；Lorentzen et al.，2002).最近，我们发现L=2附近高能质子通量的急剧下降对应着地磁SYM-H指数的下降，即使是小磁暴(-50 nT

很多卫星穿越内辐射带，都观测到MeV量级高能质子的损失.Explorer 26卫星观测到1965年4月17日磁暴期间40～110 MeV的质子通量减少(McIlwain 1966).CRRES卫星观测发现高能质子清空现象(Gussenhoven et al., 1994；Hudson et al., 1997，1998).HEO-3卫星观测发现内辐射带L≥2.3的区域，能量为27～45 MeV的质子通量在1～3天内减小(Selesnick et al., 2010).NOAA15、NOAA16和NOAA17卫星观测表明内辐射带高能质子(35～500 MeV)对大磁暴(|Dst|>200)存在响应(Zou et al., 2011).上述高能质子损失的研究主要集中在几十天以上的时间尺度，一些则是几个月的时间尺度(Lorentzen et al., 2002).最近，Van Allen Probes卫星发现内辐射带中存在与地磁活动时间尺度相类似的高能质子快速响应(Xu et al., 2019).

目前磁暴期间内辐射带高能质子的损失机制主要包括磁流体波的影响(Dragt，1961；McIlwain，1965)，磁力线曲率散射(Selesnick et al., 2010；Zou et al., 2011；Engel et al., 2015，2016)和磁场扰动(Anderson et al.,1997；Young et al.,2002；Tu et al.,2014).这些机制都是非绝热过程导致高能质子损失，然而绝热效应也能导致磁暴期间粒子通量的下降，被称为“Dst效应”.目前绝热变化主要研究外辐射带相对论电子通量的减少(Li et al., 1997；Kim and Chan, 1997).内辐射带区域高能质子的绝热通量变化关注很少.

Van Allen Probes-A观测结果表明，无论磁暴强度如何，高能质子通量的急剧下降都伴随着SYM-H指数的相应下降，并呈现一对一的对应关系(Xu et al., 2019).因此，人们很自然地会问：磁暴期间高能质子通量变化是否是磁场变化的绝热响应引起？主相期间通量下降和恢复相通量增加有多少是由于绝热效应造成的？有多少磁暴事件高能质子减小经历了完全绝热过程？本文的主要目的是定量评估完全绝热变化在多大程度上可以解释磁暴期间观测到的内辐射带高能质子(18.5～31.2 MeV)通量的急剧下降和恢复.

1 数据与方法

1.1 数据来源

范艾伦任务包括两颗卫星，Van Allen Probes-A和Van Allen Probes-B.卫星轨道的近地点距离地球约600 km、远地点距地球约30500 km.卫星在赤道面上覆盖的径向距离范围为L=1.1～6.0的区域，地磁纬度范围从-20°到20°.高能质子通量数据来自Van Allen Probes上相对论电子质子仪器(REPT)(Baker et al., 2012).本文详细介绍高能质子(能量为18.5～24.0 MeV和24.0～31.2 MeV)的变化，所有高能质子通量通过归一化方程(Xu et al., 2019中方程1和2)投射到赤道平面.

1.2 定量评估完全绝热效应

磁暴期间随着环电流的增强，内辐射带磁力线拉伸，磁场减小.为了保持磁通量Φ不变，L壳指数增加.L壳指数定义为漂移壳的径向距离.下标p、m和r分别指磁暴前、磁暴主相和恢复相阶段的所有参量.磁暴前tp时刻，高能质子在磁场强度为Bp的Lp壳层中作漂移运动，动能为Ep.磁暴前高能质子通量表示为jp=j(Ep,Lp;tp).磁暴主相期间，内辐射带磁场强度从Bp减小到Bm，L壳指数从Lp增大到Lm.此外，对于第一不变量μ守恒，质子动能从Ep降低到Em.磁暴主相高能质子通量表示为jm=j(Em,Lm;tm).恢复相期间，磁场结构恢复到近似于磁暴前的状态，强度从Bm增加到Br；质子能量从Em增加到Er，质子L壳指数从Lm减小到Lr.

Liouville定理表明，粒子的相空间密度f=j/p2(其中j指粒子通量，p指粒子动量)沿其动力学路径是恒定的(Roederer，1970).它可以表示为

f(μp,Jp=0,Φp;tp)=f(μm,Jm=0,Φm;tm),(1)

Jp=Jm=0表示赤道面，所有的观测数据都会被投影到赤道平面.磁暴主相期间动能Em和磁暴前动能Ep之间的关系可以从第一绝热不变量守恒推导出来：

(2)

pp和pm分别是磁暴前和磁暴主相期间的质子动量.相对论动量p和动能E之间的关系可表示为(pc)2=E2+2mc2E.因此，

(3)

利用第三不变量守恒，磁通Φ=∮BpdSp=∮BmdSm，Lm和Lp之间的关系可以描述为

(4)

式中k0=30115.3 nT，δB是变化的磁场，将在第1.3节详细描述.

磁暴期间高能质子通量jm与磁暴前质子通量jp的关系如下：

(5)

因此，对于绝热过程，磁暴期间高能质子通量可以根据磁暴前高能质子通量乘以磁暴前后磁场强度比值Bm/Bp给出.

1.3 磁场模型

计算完全绝热效应需要引入磁场模型：包括Hilmer-Voigt对称环形电流场模型(Kim and Chan，1997)、修正的偶极子模型(Selesnick and Kanekal，2009；Tu and Li，2011)和输入参数来自太阳风动压和行星际磁场的T96磁场模型(Tsyganenko, 1996).

图1给出了T96模型(黑线)和修正的偶极子模型(灰线)预测的L=2的磁场数据.T96磁场模型的输入参数为行星际磁场(图1a和1a′)和太阳风动压(图1b和1b′).图1c和1c′分别表示地磁平静期(2013年1月3日，左)和活跃期(2015年3月17日至18日，右)T96模型和修正的偶极子场模型的结果.磁平静期，T96模型和修正的偶极子模型计算的磁场形态几乎一致，存在恒定差约为14 nT.磁活动期间，两种模型有明显差异，尤其在2015年3月17日磁暴的主相期间扩大.恢复相期间，两个磁场模型预测的差异变小.本文使用两个磁场模型来计算磁暴主相和恢复相期间的高能质子通量.

图1 磁静(2013年1月1日)和磁暴(2015年3月17—18日)期间，T96磁场模型(黑线)和修正的偶极子场(灰线)模型计算的磁场强度

修正的偶极子模型(Selesnick and Kanekal, 2009，方程8)是指地球的偶极子场Bdip叠加一个大小等于Dst指数的均匀南向磁场，可以表示为B=Bdip+δB.δB可以表示为

(6)

通常Dst指数的时间分辨率为1 h，而本文的计算需要更高的时间分辨率.众所周知，SYM-H指数与Dst指数非常相似，时间分辨率为1 min.因此，本文使用SYM-H指数代替Dst指数.

1.4 磁静期间高能质子通量模型

为了获得磁静期间高能质子通量随能量和L壳指数的分布特征，即方程式(5)中jp=j(Ep,Lp;tp)，本文使用了Van Allen Probes-A卫星REPT仪器数据，统计平均四年(2013—2016年)的高能质子通量数据，构建了地磁平静期间(AE<200 nT，Kp<2)通量随能量和L壳指数变化的模型.

八个能通道的高能质子通量数据根据不同的L壳指数间隔，从1.1～3.0，间隔宽度为0.01，时间间隔1个月，统计平均得到每个L壳指数的通量.图2a显示了赤道面三个不同时间段第一个能量通道(E=21.25 MeV)的平均高能质子通量.黑线方形和浅灰线三角形代表2013年1月和2015年5月一个月期间的平均通量值；深灰线圆形反映的是2013—2016年四年期间通量随L壳指数变化的特征.图2b显示了三个时间间隔内，在L=2.0时八个不同能量通道的平均质子通量.

如图2a所示，两条月平均的质子通量曲线(黑线方形和浅灰线三角形)与四年平均的质子通量平均值(深灰线圆形)相差很大.高能质子平均通量随着时间的推移而增加，可能是由于被捕获太阳质子的太阳调制，Xu等(2019)文中的图5和Selesnick等(2016)文中的图2都有相似的特征.如果我们选择四年平均的质子通量作为初始输入，那么通量随时间增加的趋势被忽略.四年平均通量只是一个数值，没有反映2013—2016年期间高能质子变化的任何趋势.该数值若成为2013和2014年的初始输入，其明显大于真实的观测数据，而成为2015和2016年的初始输入时，该数值又会小于观测数据，这将导致通量预测数据不准确.因此，我们选择了每月平均的质子通量数据，保留了随时间增加趋势的月平均通量数据.我们发现将每个月的磁静期间的平均通量作为磁暴前质子初始通量的输入，比四年的平均通量作为输入要准确得多.

图2 (a)三个不同时间段(2013年1月、2015年5月和2013—2016年)磁静时间(AE<200 nT和Kp<2)内辐射带高能质子通量随L壳指数的变化;(b)三个时间间隔内L=2的区域，高能质子通量随能量的变化曲线

基于磁场模型和磁静期间高能质子通量的分布，我们定量评估了磁暴前后高能质子的绝热变化.计算方法如下所述：(1)将观测到的高能质子通量归一化到磁赤道位置(Xu et al., 2019，方程1和2)；(2)根据已知参数Lm=2.0和Em=21.25 MeV，通过第一和第三个不变量守恒，求出Lp和Ep；(3)根据磁静期间高能质子通量分布，通过线性插值，追溯出相应的j(Ep,Lp;tp)；(4)根据方程5计算j(Em,Lm;tm).(5)以上四个步骤分别在修正的偶极子场和T96磁场模型中进行计算.

2 观测结果

2.1 事例研究：2015年3月17日和2016年1月20日的两次磁暴事件

为了研究磁暴期间高能质子的完全绝热效应，本文选择2015年3月17日和2016年1月20日两次磁暴事件作为典型例子.

磁暴主相定义为磁暴开始后从SYM-H指数低于-15 nT时开始，到SYM-H指数最小值的时间段.恢复阶段包括从SYM-H指数最小值开始到SYM-H指数恢复到最小值的75%的时间段.时间间隔分别由三条垂直虚线突出显示，如图3和图4所示.

图3显示2015年3月17日事件的计算结果.左侧(Em=21.25 MeV，Lm=2.0)和右侧(Em=27.6 MeV，Lm=2.0)分别表示两个不同的能量通道.值得注意的是，所有参数均由两个磁场模型计算得出.灰线代表修正的偶极子场的结果，黑线代表T96磁场模型的结果.图3第一行是根据参数Lm=2.0推导出来的磁静期间Lp壳指数.第二行是根据参数Em=21.25 MeV和Em=27.6 MeV反推得到的磁静期间高能质子的能量Ep.图3c和3c′表示磁静期间初始的高能质子通量分布j(Ep,Lp)，其中Ep和Lp显示在图3a、3a′、3b和3b′中.图3d和3d′显示两个能量段Van Allen Probes两颗卫星(灰色圆点表示Van Allen Probes-A和灰色*表示Van Allen Probes-B)的观测数据，以及利用两个磁场模型理论计算的完全绝热通量(灰线和黑线).图3e和3e′显示两个不同磁场模型在L=2.0时的磁场强度，从图3e和3e′可以看出，T96模型的磁场强度变化大于修正的偶极子模型结果.图3f表示地磁活动SYM-H指数.图4描述了2016年1月20日的磁暴事件，格式与图3完全相同.

图3 2015年3月17日磁暴期间，根据两个磁场模型计算在Lm=2.0时两个能量段(Em=21.25 MeV，左列；Em=27.6 MeV，右列)的高能质子绝热效应的通量变化

图4 2016年1月20日磁暴期间理论计算绝热效应导致的通量变化(jm)和观测结果.格式与图3相同

观测得到的高能质子通量数据在每个轨道周期内进行了统计平均处理.理论计算和观测结果都发现磁暴期间高能质子通量随SYM-H指数变化而变化，而且具有相同的时间尺度.对于2015年3月17日磁暴事件，主相阶段，两种磁场模型计算出的完全绝热通量下降与范艾伦卫星的观测结果相当吻合.而磁暴前和恢复相阶段，T96模型计算得到的结果比修正的偶极子模型的结果要好得多.对于2016年1月20日磁暴事件，T96模型得到的结果在整个磁暴时间段都比修正的偶极子模型好得多.修正的偶极子模型得到的理论预测数据略大于T96模型的预测数据.我们使用均方根预测误差(RMSPE)估计观测数据和理论预测数据之间的误差大小:

yobsandypre分别代表观测数据和理论计算结果.本文计算了两个能量通道和两个磁暴不同发展阶段(磁暴主相和恢复相).均方根误差都小于1%(主相期间两个能量通道分别为0.52%和0.62%；恢复相阶段两个能量通道分别为0.82%和0.9%).

2.2 统计分析：2013—2016年期间67次磁暴事件

利用磁场模型和范艾伦卫星数据建立的磁静期间高能质子通量作为输入，可以计算Xu等(2019)表1中67次磁暴期间完全绝热效应所导致的高能质子通量变化.图5显示了两颗范艾伦卫星(分别用黑色圆圈和星号表示)在磁暴主相和恢复相阶段高能质子通量的观测数据以及基于修正的偶极子(左侧)和T96(右侧)模型理论计算得到的结果.上面两排(图5a、5a′，5b和5b′)属于E=21.25 MeV能量通道，下面两排(图5c、5c′，5d和5d′)属于E=27.6 MeV能量通道.另外磁暴主相(图5a、5a′、5c和5c′)和恢复相(图5b、5b′、5d和5d′)阶段进行了分开的统计调查.

如图5所示，横坐标是观测数据，纵坐标是利用方程(3)理论计算的高能质子通量.黑线表示观测值和预测值的拟合关系Y=aX，原则上观测值与理论结果的关系为Y=X，实际上根据黑线可以看出有偏差.对于修正的偶极子模型，主相期间两个能量段的参数a分别为0.932和0.951，恢复相期间两个能量段的参数a分别为0.961和0.964.预测的通量值与观测值之间的相关系数都大于0.87.对于T96模型，两个能量段的参数a在主相期间为0.888和0.912，在恢复相阶段分别为0.900和0.950.预测通量与观测通量的相关系数分别大于0.84.从图5b、5b′、5d和5d′可以看出，在恢复相阶段，修正偶极子模型的预测通量比T96模型的预测通量分布更加分散.正如我们所预期的，两个阶段预测的通量与观测值高度相关.因此，在磁暴主相和恢复相期间，完全绝热效应可能是内辐射带高能质子减少和恢复的主要原因.

从图5也可以看出，卫星观测通量数据与理论预测值之间存在差异.完全绝热效应贡献了大约90%的高能质子通量下降和恢复，这意味着可能还涉及一些非绝热损失机制，如上所述的低频电磁波动和磁力线曲率散射(McIlwain, 1965; Anderson et al., 1997; Young et al.,2002; Tu et al.,2014; Engel et al., 2015，2016).在磁暴恢复阶段，非绝热效应对高能质子通量的恢复可能起着重要作用，不应被忽视.因此本文进一步分析了67次磁暴期间相空间密度的变化，试图找出有多少磁暴经历了绝热过程，有多少磁暴不能用绝热效应来描述.本文使用的相空间密度数据直接从Van Allen Probes数据网站下载：https:∥rbspgway.jhuapl.edu/psd.

2.3 2013—2016年期间磁暴事件相空间密度分析

图6表示磁暴三个不同阶段的高能质子通量和相空间密度随L*壳指数的分布.深灰色线表示磁暴前，黑线表示主相，浅灰色线表示恢复相.左右两列分别表示两个能量段(E=21.25 MeV和E=27.6 MeV)从上至下分别表示高能质子通量(图6a和6a′)、K=0.08 G1/2RE和K=0.11 G1/2RE时的相空间密度(图6b、6c和6b′、6c′).图6d和6d′表示在磁暴前和磁暴后的高能质子通量.如图6a和6a′所示，磁暴主相期间，高能质子通量在L壳指数1.75和2.5之间显著下降.磁暴恢复相阶段，质子通量恢复到磁暴前的水平.深灰色线和浅灰色线几乎一致，这说明粒子经历了绝热过程.图7表示发生于2016年1月20日的磁暴，格式与图6相同.从图7发现2016年1月20日磁暴期间发生高能质子通量的变化属于绝热效应.

图6 2015年3月17日磁暴三个时间段的通量分布(a和a′)，K=0.08 G1/2RE(b和b′)和K=0.11 G1/2RE(c和c′)时的相空间密度，以及磁暴前后高能质子通量的变化率(d和d′)

图7 2016年1月20日磁暴事件的通量、相空间密度以及磁暴前后质子通量的变化，格式与图6相同

图8显示2015年12月20日磁暴期间高能质子的非绝热变化.图8的格式与图6和图7相同.磁暴主相期间，几乎整个内辐射带的高能质子通量急剧下降.五天后，SYM-H指数恢复到磁暴前水平；L<1.75的内辐射带区域，质子通量几乎恢复到风暴前的水平；然而L壳指数在1.75～2.5之间的区域，质子通量恢复到磁暴前水平的60%左右.相空间密度的变化也表现出与观测通量相似的特征.

图8 2015年12月20日磁暴事件的通量、相空间密度以及磁暴前后质子通量的变化，格式与图6和图7相同

磁暴前后高能质子通量的统计结果如图9所示.根据磁暴前(图6a深灰色线)后(图6a浅灰色线)高能质子随L*变化的剖面，定义两者的变化率(如图6d和图6d′)来讨论磁暴前后高能质子变化究竟是不是由绝热效应造成的.图9所示67次地磁暴中，两个能量段E=21.25 MeV和E=27.6 MeV的高能质子有56次事件经历了绝热变化，a>0.95.因此，本文认为在大多数(56/67)地磁暴期间，绝热变化在内辐射带是比较常见的.

图9 67例磁暴事件前后高能质子通量变化率的统计结果

如图10所示，(a、b、c)分别表示SYM-H指数和两个质子能量通道(21.25 MeV和27.6 MeV)，与Xu等(2019)的图3相似.图10d和10e是两个不同L*(黑色代表L*=2.0，灰色代表L*=2.3)在赤道面(Jm=0)u=535 MeV/G和u=700 MeV/G的相空间密度，它们对应的能量为21.25 MeV和27.6 MeV.垂直虚线表示2013—2016年间的67次磁暴事件.从图10d和10e可以看出，对于中小磁暴，磁暴前后的相空间密度变化较小，说明绝热过程控制着磁暴前后的通量变化.对于某些磁暴，如图8所示磁暴前后相空间密度变化很大，意味着强烈的地磁活动可能涉及非绝热过程，这类高能质子通量减小或清空的物理机制有待进一步研究.

(图10续)

3 总结

本文发现内辐射带外区能量小于40 MeV的高能质子通量在磁暴主相期间显著减少，并随着SYM-H指数的恢复而恢复.高能质子通量的变化与磁场和SYM-H指数的变化具有相似的时间尺度.基于这些观测结果，利用刘维尔定理和第一绝热不变量和第三绝热不变量守恒，检验了磁暴主相和恢复相阶段内辐射带高能质子的完全绝热效应.利用事例分析和统计研究两种方式对绝热效应进行了定量评估.一种是利用初始通量来追踪磁暴期间的质子通量，另一种是分析磁暴前和磁暴后的通量变化和相空间密度.2015年3月17日和2016年1月20日两个事例研究发现，计算出的完全绝热通量下降与Van Allen Probes-A卫星在主相和恢复阶段的观测结果相当吻合.相比于修正的偶极子场模型，T96磁场模型的结果更吻合主相和恢复相阶段的观测结果.理论计算和观测到的磁暴期间高能质子通量随SYM-H指数变化而变化的时间尺度是一致的.

本文计算了2013—2016年发生的67次磁暴期间的高能质子通量.磁暴主相和恢复相阶段，两个磁场模型的相关系数分别大于0.84.高相关系数表明，完全绝热效应是内辐射带高能质子磁暴主相期间减少，在恢复相期间增加的主要助力.相空间密度分析表明大约83%(56/67)的磁暴前后，高能质子的相空间密度始终保持不变，与利用磁场模型预测高能质子通量的事例研究和统计结果一致.因此大部分(56/67)磁暴事件中完全绝热效应对磁暴期间高能质子通量下降和恢复起着主要贡献.内辐射带绝热行为的原因是在较低L壳指数的区域有较强的磁场，需要更强烈的外部扰动，以及高质量质子(与电子相比)对振荡的响应较慢.

此外，我们发现大于40 MeV的五个能量段的高能质子在1.75

致谢感谢范艾伦卫星科研团队提供数据.范艾伦卫星高能质子和空间相密度数据来源于(ftp:∥cdaweb.gsfc.nasa.gov/pub/data/rbsp/和https:∥rbspgway.jhuapl.edu/psd).感谢京都世界地磁数据中心提供SYM-H指数数据(http:∥wdc.kugi.kyoto-u.ac.jp/).I.R.感谢NASA经费NNN06AA01C的支持.