EP和α “双追踪” GRB 131231A的能谱特性研究

王 晨，彭朝阳，陈家明

(云南师范大学物理与电子信息学院，云南 昆明 650500)

人类研究伽马暴已经有50多年的历史，伽马暴的瞬时辐射在几个基本问题上依然让人感到困惑，一些非常基础的问题是喷流成分(重子主导的火球或是坡印廷流量(Poynting-flux)主导的流出物)，它决定或至少强烈影响能量耗散机制(激波或磁重联)、粒子加速机制(热主导或磁主导)，特别是辐射机制(同步辐射)或光球准热辐射的康普顿化[1]。伽马暴辐射有瞬时辐射和余辉两个主要阶段。瞬时辐射阶段通常持续几秒到几分钟。伽马暴瞬时辐射具有丰富的观测特征，但是其中的大多数没有合理的解释。伽马暴瞬时辐射阶段的光谱成分主要来自两种起源：同步辐射(非热成分)和光球层(热成分)[2]。过去的研究表明，伽马暴能谱中可能存在3种基础成分[3]：(1)能用Band函数拟合的非热成分；(2)能用普朗克函数拟合的热成分；(3)出现在能量较高区域的幂律成分。目前已观测到的伽马暴中能发现两种或两种以上的成分[4]。热成分主导的伽马暴非常稀少，所以对于暴的能谱处理时多采用组合模型，常用的基础模型是Band函数或者截止幂律(Power law with exponential cut off; CPL)函数。多数情况下，组合模型能得到更优的参数值，有助于我们推断伽马暴能源机制和中心引擎。在BATSE探测器时期，Ep显示出几种不同的演化：(1)“从硬到软” 的趋势，无论流量如何升降，它都是单调递减的[5-6]；(2)“流量跟踪” 趋势，即参数的演化跟随流量上升或下降[7]；(3)“软到硬” 的趋势或混乱的演变[8]。2008年费米(Fermi)卫星升空后，大量的伽马暴观测数据证实前两种模式占主导地位： “由硬到软” 约占三分之二， “流量跟踪” 约占三分之一[9]。这些演化模式的物理起源仍然没有解决，尽管文献中已经提出了一些方案[10-11]。与Ep演化相比，α演化更加混乱，因此研究和物理解释相对较少。文[9]发现GRB 131231A的Ep和α演化都表现出 “通量跟踪” 行为，并且定义为光谱演化的 “双跟踪” 模式[9]，但文[9]并没有加入热成分进行详细的能谱分析。既然很多研究表明伽马暴能谱中包含热成分，为了更详细地研究光谱参数EP和α“双追踪” 的能谱特性，我们使用贝叶斯分区并使用蒙特卡罗模拟方法处理GRB 131231A。在伽马暴的研究中，研究人员发现次主导的热成分，一般出现在Band谱[12]的左侧[13]，即热成分一般出现在低能Band谱区域。针对这种现象，文[9]在探究8个具有较明显热成分的伽马暴中热成分对于非热成分几个特征参数值的影响时，都选用截止幂律函数作为基础模型与普朗克函数(Blackbody, BB)模型组合进行探究[2]。对于绝大多数伽马暴的时间分辨谱，截止幂律函数模型是适合的[14-15]。然而在一些例子中，Band模型的高能幂律指数β显著提高拟合。在同步辐射的情况下，高能幂律指数β提供了有关粒子加速和激波的本质，而在光球辐射的情况下，它提供了有关辐射区域能量耗散的信息[16]。考虑到这种情况，我们以统计值BIC(Bayesian Information Criterion)作为判断标准[17]，对伽马暴的每个拟合区间选用最优模型的数据，并研究使用最优模型之后GRB 131231A能谱参数的变化情况。

1 数据处理和分析方法

GRB 131231A由费米卫星中的伽马暴探测器(8 keV～40 MeV)于2013的12月31触发。本文中，我们使用贝叶斯分析包，即最大似然工作框架(3ML)分析所有伽马暴的光谱。该工具包采用马尔可夫链蒙特卡罗理论(Markov Chain Monte Carlo, MCMC)对时间分辨谱进行处理,这种方法的优点详见文[1]。我们从网站(https://fermi.gsfc.nasa.gov/)获取伽马暴探测器[18]三种类型的数据，CTIME文件、CSPEC文件和TTE文件，选择具有最小时间分辨率(2 μs)的TTE文件进行数据分析。在贝叶斯统计的马尔可夫链蒙特卡罗拟合例子中，为了确定能谱拟合的最优模型，我们选取BIC[17]作为判断标准，通过比较不同模型BIC的差值(ΔBIC)确定最佳模型。一般来说，Band模型[12]和截止幂律函数模型[19]在伽马暴拟合时是两种较为优越的模型，我们把这两种模型选作基础模型，分别与普朗克函数模型组合，在伽马暴时间分辨谱的每一个区间经过模型对比，进而选出适应该区间的最优模型。为了更好地研究光谱随时间的演化，对光谱进行时间分辨谱分区是一种有效的研究方法。文[20]研究了4种对光变曲线分区时间片的工具，其中信噪比(S/N)和贝叶斯分区(BBs)[21]是两种最有效的分区方法。这两种方法各有优缺点，传统的信噪比分区可以确保每个时间段能分到足够多的光子做光谱拟合，但是有时它会破坏物理结构。关于贝叶斯分区的优缺点详见文[20]，我们采用贝叶斯分区，但是依然使用信噪比作为一种筛选方法剔除所得区间中不满足某一标准的区间(本文选取信噪比S/N≥20的区间，因为信噪比S/N<20的区间拟合时出现较大误差[1])，共同使用这两种方法，我们可以得到更有研究意义的区间。表1给出了该暴的相关参数，注：CPL1代表ΔBIC≤-8的时间分辨谱区间使用截止幂律函数模型，CPL2代表-8<ΔBIC<8的时间分辨谱区间使用截止幂律函数模型。图1给出了一个分辨谱的拟合结果。

表1 伽马暴参数及使用的探测器编号Table 1 Gamma-ray burst parameters and detector numbers used

图1 一个黑体成分最显著的计数率时间分辨谱(时间区间为24.59～25.86)拟合结果。

2 分析结果

我们使用Band函数模型、康普顿函数模型和普朗克函数模型对GRB 131231A进行拟合，得到它的时间分辨谱数据。从时间分辨谱的ΔBIC来看，有些暴的某些时间段中截止幂律函数模型的拟合结果优于Band函数模型，有些结果相反。综合考虑之后我们把每个暴的每个时间分辨谱区间进行最优模型选择，选择方法是用每个区间单独用截止幂律函数模型拟合时产生的BIC减去单独用Band函数模型拟合时产生的ΔBIC作为筛选标准，当ΔBIC≤-8时我们认为该区间更适合使用截止幂律函数模型，当ΔBIC≥8时我们认为该区间更适合使用Band函数模型，当-8<ΔBIC<8时我们认为该区间更适合使用截止幂律函数模型(当两种模型中任何一种模型优势不明显时，一般选用参数少的模型)。然后再使用该区间选定基础模型的BIC与基础模型加普朗克函数模型的组合模型的BIC差值是否≥8判断该区间是否有热成分，最后筛选信噪比S/N≥20的区间作为研究对象。

光谱分析得到单个时间区间的时间跨度远小于整体的时间跨度，我们假设在每个时间区间中参数的变化是准静态的，用时间区间中值时刻参量的特征代表整个时间区间参数的行为[22]。

2.1 参数的分布

图2显示了Best模型参数分布，包括α,Ep和Fp。表2列出了分布的平均值和标准偏差。我们对各个参数分布采用高斯函数核进行核密度估计，结果如表2。未加入黑体成分和加入黑体成分时，低能光子指数α分布的均值分别为-1.06 ± 0.29和-1.00 ± 0.29，加入黑体成分后，α轻微变硬，但与典型值α～-1基本一致，且没有超过所谓的同步 “死线”(-3/2<α<-2/3)[23]。在加入黑体成分后，峰值能量EP的核密度曲线峰值向右偏移，均值为199.84 ± 161.61，与未加入黑体成分相比，峰值能量Ep变硬。加入黑体前后Fp没有明显变化，平均值仅有微弱变化，这些结果和前人的研究结果一致。

图2 Best模型的参数α，EP和FP分布Fig.2 The distribution of the parameter values α, EP and FP

表2 模型参数均值与标准偏差

2.2 加入黑体前光谱参数的演化

我们用本文分析的数据重新分析未加入黑体成分的峰值能量和低能谱指数的演化情况。图3展示了贝叶斯区间中S/N≥20的时间分辨谱区间对应的数据点。GRB 131231A的参数EP和α表现为较好的 “流量跟踪” 趋势，本文的结果与文[9]的结果一致。

图3 GRB 131231A最优模型的时间分辨谱参数EP和α的分布Fig.3 The distribution of the parameter values EP and α of the time-resolved spectrum of the optimal model of GRB 131231A

2.3 加入热成分后光谱参数的演化

加入黑体成分后，我们找出探测到热成分的16个时间分辨谱区间发现，Ep和α的演化趋势和未加入热成分的趋势一致，依然表现出双追踪行为。这两个参数发生了变化，α和Ep轻微变硬，但从图4中我们能发现这两个参数变化的细节，GRB 131231A中有热成分的区间EP值的整体变化并不一致，如图4(a)。我们发现GRB 131231A中EP在峰值之前，加入普朗克成分后变硬(值变大)；峰值之后，加入普朗克成分后变软(值变小)。其中EP变硬的点有6个，均分布在峰值之前；EP变软的点有10个，均分布在峰值之后。GRB 131231A中α的整体变化也不一致，如图4(b)。有一个点分布在同步辐射死线之上，其余点分布在同步辐射死线之下，加入普朗克成分后，α出现了明显的变化，在峰值之前，加入普朗克成分后α值变软；峰值之后，加入普朗克成分后变硬。在加入黑体成分前，GRB 131231A中仅有最后两个点的α值小于-3/2；加入热成分之后，脉冲峰值处有两个点的α值上升到超过同步辐射死线，最后两个点的α值略有上升，但是依然小于-3/2。

图4 (a)GRB 131231A中探测到热成分的时间分辨谱区间中EP和α的分布；(b)上下两条红色虚线分别代表

2.4 光谱参数的相关性

相关分析在理解伽马暴物理学中起着重要作用，因为它提供了一些线索帮助我们揭示伽马暴的本质。我们探究了参数logEP-logFP，logFP-α和logEP-α的关系，表3给出了具体相关系数。相关强度由斯皮尔曼相关系数r给出，强相关性r>0.7，弱相关性r<0.4。文[7, 24-25]揭示了logF-α的关系主要有3种类型：非单调关系(包含正和负幂律段，在峰值流量处具有明显的间断)、单调关系(由单个幂律描述)和没有明显的趋势。GRB 131231A 3种模型中logEP-logFP的相关性，Best为r=0.426,p=0.037 8；Band为r=0.687,p<10-4；CPL为r=0.566,p<10-4，最好的是Band模型。3种模型中logEP-α的相关性，Best为r=0.412,p=0.045 3；Band为r=0.70,p=0.000 14；CPL为r=0.731,p<10-4，最好的是Band模型。3种模型中logFP-α的相关性，Best为r=0.765,p<10-4；Band为r=0.655,p=0.000 52；CPL为r=0.756,p<10-4，最好的是Best模型。在Band和CPL模型中logEP-logFP，logFP-α和logEP-α都表现出极强的相关性。当我们使用Best模型时，logEP-logFP和logEP-α的相关性均变弱；logFP-α相关性强于Band函数和截止幂律函数。

表3 GRB 131231A在不同模型中EP-Fp，EP-α和Fp-α参数关系(斯皮尔曼(Spearman)相关性系数r和置信度p)

2.5 热成分分析

为了检查能谱中的热成分，我们先寻找对应于时间分辨谱的最优模型，用最优模型加入普朗克函数成分拟合，查看加入黑体前后统计量的差值ΔBIC，如果ΔBIC<-8说明光谱中显著含有热成分。图5(a)给出了ΔBIC在脉冲中的演化和分布图，图5(b)给出了热成分流量与总能量流量的比率。

GRB 131231A具有明显的脉冲上升阶段和下降阶段。我们对GRB 131231A中出现热成分的区域进行统计，从图5(a)统计的结果来看，ΔBIC<-8的光谱主要出现在对应高光子计数的光变曲线的位置，说明热成分主要存在于脉冲前期，并且随着时间逐渐减少。图5(b)表明，热成分流量与总能量流量的比率多数分布在0.1以下，比率较高的区域主要分布在高脉冲区域。

图5 (a)ΔBIC随脉冲的演化和分布图，绿色和红色虚线分别表示ΔBIC等于0和-8的位置；(b)探测到热成分的光谱中热成分流量与总能量流量的比率，红色虚线代表比率为1

2.6 与热成分相关的参数

现有的研究表明，热成分起源于光球层，根据高纬辐射解释，最高温度和最强热流最初从光球层辐射轴发出，并朝向观测者。因此，原则上来说，光子辐射位置的半径可以确定[26]。按照文[26]的理论，对于一个已知红移的伽马暴，知道该暴的观测温度和热能量流量时，我们可以计算流出物的初始半径r0(从中心引擎到喷流喷射处的距离)，对于rph>rs的情况可以更进一步计算伽马暴的饱和层半径rs、光球层半径rph和流出物滑行运动阶段洛伦兹因子Γ。根据文[22]中(1)式可以判断光球层和饱和层的相对位置，我们发现GRB 131231A中有热成分区间均满足(FBB/FNT)-3/2>1，判定这个暴的所有壳层的光球层位于饱和层之上(rph>rs)，应用文[26]中的理论计算这个暴的初始半径r0、饱和层半径rs、光球层半径rph和流出物滑行运动阶段洛伦兹因子Γ。

我们能得到与热成分相关的参数R，R可以解释为辐射区的有效横向尺度[22]，R为常数意味着光球的有效辐射区域与时间无关。我们可以通过

(1)

计算得到伽马暴的R值，其中，σSB为斯特藩-玻尔兹曼常量，单位为erg·s-1·cm-2·K-4；FBB为观测到的黑体能量流量，单位为erg·s-1·cm-2；T为黑体的绝对温度，单位为keV；z为伽马暴的宇宙学红移；dL为光度距离，单位为cm；Γ为流出物滑行阶段的洛伦兹因子。

如图6(a)，在脉冲起始阶段R表现为快速上升，在22 s后R没有表现出单调递增的趋势，变化较为平缓。流出物速度方向与视线的夹角为θ，在脉冲峰值时刻之后，θ=0区域的辐射减弱，高纬辐射起主导作用[26]，此时光球的有效辐射区域与时间无关。

为了分析温度随时间的演化，文[27]依据文[28]的拐折幂律模型派生了一种模拟温度的拐折幂律模型。拐折幂律的形式如文[27]中(6)式，我们使用它拟合该暴探测到热成分区间的温度时间，如图6(b)。我们得到对应于(6)式的参量拟合值分别为Φ=2.77 ± 2.44，t0=24.59 ± 0.95，ξ=-4.88 ± 2.49和δ=0.0028 ± 1.08；由拟合参量我们得到拐折幂律函数的上升和下降阶段幂律指数分别为a=7.65 ± 3.48和b=-2.10 ± 3.48。

图6 (a)R 随脉冲的演化图；(b)探测到热成分的光谱中温度随脉冲的演化，蓝色虚线为拟合线

GRB 131231A的洛伦兹因子呈现先上升后随时间衰减的演化趋势，起始值为157 ± 0.5，峰值为385.6 ± 0.34。它能用拐折幂律较好地拟合，如图7(a)，上升阶段的幂律指数α1=1.37 ± 0.39，下降阶段的幂律指数β1=-1.38 ± 0.39。文[22, 29]发现在具有显著光球成分的伽马暴中，洛伦兹因子随时间单调减小，Γ∝(F/R)1/4Y1/4。在脉冲上升阶段，总流量F和R的增加速率几乎相同，洛伦兹因子接近固定值或仅随时间微弱衰减。我们得到的洛伦兹因子呈现为先增加后衰减，与传统光球模型预言的洛伦兹因子在初始阶段近似为常数不符。

图7 GRB 131231A的Γ，r0，rs和rph的演化Fig.7 The evolution of Γ, r0, rs and rph of GRB 131231A

初始半径表示喷流开始加速。在GRB 131231A中r0从初始值(3.6 × 106± 1.56 × 106)Y3/2cm开始随光变曲线上升，在脉冲峰值处(22 s)达到最大值(1.39 × 108± 5.17 × 107)Y3/2cm，但是并没有很好地追踪光变曲线；脉冲峰值后，r0又呈现上升的趋势，这可能是由于中心引擎的再活动导致；在脉冲后期，r0下降至(6.9 × 105± 1.02 × 106)Y3/2cm，低于初始值。r0的变化量级分布在106～108，这个演化范围在黑洞的视界半径与Wolf-Rayet前身星的核心大小之间，表明r0与喷流和前身星之间的相互作用有关[22,29]。

饱和层半径rs的演化与初始半径在脉冲中的行为相似，在数值上相差近3个数量级。rs初始值为(5.7 × 108± 2.46 × 108)Y-5/4cm，最大值为(4.72 × 1010± 1.76 × 1010)Y-5/4cm，演化过程中的最小值为(1.2 × 108± 1.83 × 108)Y-5/4cm。

3 结论与讨论

GRB 131231A是一个比较明亮且较为接近单脉冲的长暴。我们应用贝叶斯方法，利用经验光谱模型和物理模型，即截止幂律函数模型、Band函数模型和普朗克函数模型，进行详细的时间分辨谱分析，并在每个时间分辨谱区间上进行模型比较，确定适合某一区间的最优模型，我们总共获得了24个光谱并研究了它们的性质，这些光谱的最佳模型的拟合参数见表4。这个伽马暴的EP和α都表现为较好的追踪行为。EP的流量追踪行为可能与两种同步辐射模型有联系[8,30]。图2(b)中显示α均分布于同步辐射死线之下，整体分布非常符合同步辐射模型。文[9]指出，在上升阶段，α变得更硬，表明电子正在从快速冷却区向慢速冷却区演化[9]。在衰减阶段，α变得更软，宽脉冲的衰减阶段可能由 “曲率效应” 控制[31-33]。

表4 伽马暴131231A最优模型拟合的时间分辨谱Table 4 The best fitting results of models to fit the time-resolved spectrum of GRB 131231A

续表4

通过最优模型与普朗克函数模型的共同拟合，我们探测到16个有明显热成分的时间分辨谱区间，该暴不是热主导的伽马暴。文[9]系统研究了8个含有热成分伽马暴的能谱特性[1]，对于非热成分主导的伽马暴，加入热成分后，α值变软，Ec变硬，能量流量Fp在不同模型中没有明显的不同。对于GRB 131231A，在光变曲线峰值之前，加入普朗克成分后EP变硬，峰值之后，加入普朗克成分后EP变软；在光变曲线峰值之前，加入普朗克成分后α变软，峰值之后，加入普朗克成分后变硬。GRB 131231A具有特殊的光谱演化。

伽马暴喷流的中央引擎很可能是一个混合系统，既有一个热的火球成分，又有一个冷的坡印廷流量成分[34]。一些明亮的暴(例如GRB 090902B)在脉冲开始处或峰值附近表现出类似热的时间分辨谱，在随后的脉冲中出现非热的时间分辨谱，这可能表明喷流成分在逐渐变化[35]。文[34]在解释伽马暴喷流成分的起源方面做了更全面的描述，即混合喷流模型。文[9]使用此模型分析了8个费米伽马暴[1]，发现大多数爆发可以用混合喷流模型很好地解释。GRB 131231A的低能光子指数α的值不超过所谓的同步辐射死线；但加入热成分后，拟合统计量明显改善，其中ΔBIC<-8的光谱主要集中在脉冲开始阶段和峰值阶段(10～26 s)，随后热成分逐渐减少，说明伽马暴可能从物质主导火球过渡到磁主导的坡印廷流量；在脉冲后期又出现热成分，这可能由中心引擎再活动导致。

虽然温度也能用拐折幂律拟合，但是与文[27]拟合伽马暴中热脉冲得到的上升阶段和下降阶段的指数值均存在较大差异。

我们计算了与热成分相关的参数Γ，r0，rs和rph，得到的洛伦兹因子先增加后衰减，与传统光球模型预言的洛伦兹因子在初始阶段近似常数不符[22,29]，之后以指数衰减，这可能是因为内激波模型导致[36]。另一种解释是孤立磁星诞生的情况，此时洛伦兹因子随时间增加[37]，但是GRB 131231A的光度远大于磁星模型的光度，所以磁星模型的设想不太可能成立。我们得到r0，rs和rph的变化范围与具有显著光球成分的伽马暴GRB 081224A，GRB 090719A和GRB 100707A等[22,29]变化范围吻合。其中，r0的演化范围在黑洞的视界半径与Wolf-Rayet前身星的核心大小之间，表明r0与喷流和前身星之间的相互作用有关[29]。rs的演化行为与r0近似，在数值上相差3个数量级以内。在脉冲起始至达到峰值阶段(22 s)，光球层半径rph呈现跟随脉冲上升(追踪)；脉冲峰值后，rph快速下降，变化范围在一个量级内，之后仅微弱变化。脉冲的上升阶段，rph随脉冲上升，Γ也随时间增加，这要求L0或Μ以更快的速率增加，而重子载量的增加，导致与它相关的电子增加了不透明度，从而增加了光球层半径。在脉冲衰减阶段，Γ和L0都减小，而rph保持常数或仅微弱变化，这种状况下，L0∝Γ3(或Μ∝Γ2)[29]。

GRB 131231A的能谱特征行为表明这个暴有一些特殊性，需要寻求合适的理论解释，值得我们进一步深入研究。