引力波引起的脉冲星计时残差模拟与分析

童明雷，丁勇恒，赵成仕，高峰，高玉平

（1.中国科学院国家授时中心，西安710600；2.中国科学院时间频率基准重点实验室，西安710600；3.中国科学院星系宇宙学重点实验室，合肥230026；4.中国科学院大学，北京100049）

引力波引起的脉冲星计时残差模拟与分析

童明雷1，2，3，丁勇恒1，4，赵成仕1，2，高峰1，4，高玉平1，2

（1.中国科学院国家授时中心，西安710600；2.中国科学院时间频率基准重点实验室，西安710600；3.中国科学院星系宇宙学重点实验室，合肥230026；4.中国科学院大学，北京100049）

为了能将引力波引起的脉冲星计时残差从总的脉冲星计时残差中分离出来，我们介绍了3种形式的引力波（单一持续的引力波、单一瞬时的引力波和随机引力波背景）及其对脉冲星计时残差的影响，并利用脉冲星计时处理软件Tempo2进行了模拟，对模拟得到的这三种引力波源引起的脉冲星计时残差进行了分析比较。该工作在脉冲星时间尺度的建立及其应用方面具有参考意义。

脉冲星；引力波；计时残差

0 引言

毫秒脉冲星（自转周期通常小于20 ms）的高度自转稳定性使之成为自然界天然的时钟，其周期导数一般为10-20。以脉冲星自转建立的时间尺度称为脉冲星时，而由多颗自转稳定、计时噪声小的毫秒脉冲星通过适当的算法可建立综合脉冲星时间标准，这是脉冲星计时的一个重要应用。所谓脉冲星计时就是确定脉冲到达时间（TOA）的过程。由于单脉冲的强度不稳定，不适合用来计时，而将几百到几千个单脉冲作周期折叠后平均，可得到稳定的积分脉冲轮廓。将观测得到的积分脉冲轮廓与标准脉冲轮廓作相关处理，可以确定脉冲到达时间（TOA），由观测者的参考原子钟来记录，然后经过一系列的钟差修正转换到地球时（TT），再经过涉及Roemer延迟、Shapiro延迟、Einstein延迟和大气延迟等的计时模型转换到太阳系质心坐标时（TCB）。然后将其与脉冲星钟模型在太阳系质心（SSB）处预报的TCB作差，得到脉冲星计时残差。脉冲星钟模型用相位可表示为：引起脉冲星计时残差的原因有很多，比如TOA的测量误差、脉冲星内部的物理过程引起的计时噪声、星际介质的色散、原子钟钟差、宇宙中的引力波、地球历表误差及计时模型误差等。多颗脉冲星组成计时阵（PTA），目前国际上有澳大利亚的PPTA、欧洲的EPTA以及北美的NANOGrav，这三者联合起来组成了国际脉冲星计时阵IPTA。PTA使脉冲星计时得到广泛的应用，例如，可以探测引力波[1]、测量太阳系行星质量[2]、建立脉冲星时间标准[3]。未来通过太空X射线望远镜对多颗脉冲星的观测，还可以应用到深空自主导航领域。

宇宙中的引力波经过脉冲星和地球时会产生额外的脉冲星计时残差。若以探测引力波为目的，引力波就是信号，可以通过多颗毫秒脉冲星计时残差的相关性提取引力波信号。虽然原子钟钟差和地球历表误差也会引起脉冲星计时残差的相关性，但是与引力波引起的脉冲星计时残差的相关性还是有区别的。原子钟钟差对观测的所有脉冲星计时残差影响相同，呈单极性；对于位于空间相同方向的两颗脉冲星，地球历表误差对其计时残差影响符号相同，对位于空间相反方向的两颗脉冲星，地球历表误差对其计时残差影响符号相反，大小则与脉冲星的黄纬有关，呈现偶极性；引力波对相对于测站张角为0°或180°分布的两颗脉冲星的计时残差影响符号相同，而对张角为90°的两颗脉冲星的计时残差影响符号相反，呈四极性。但对于以建立脉冲星时间标准为目的，引力波就成了噪声，很好地认识这种噪声进而找到消除或减弱的方法，便可以提高脉冲时时间标准的精度。因此，无论出于何种目的，分析各种形式的引力波如何影响脉冲星计时残差都是有意义的。本文主要讨论3种形式的引力波：单一持续的引力波、单一瞬时的引力波和随机引力波背景。利用国际上被广泛应用的脉冲星计时处理软件Tempo2[4-6]，模拟各种形式的引力波，得到不同的计时残差。分析这些计时残差的特征，有利于将引力波的影响纳入到脉冲星计时模型中去，建立更高精度的脉冲星时间标准。Tempo2是较为精确的脉冲星计时处理软件，从TT转换到TCB精度优于1 ns。关于Tempo2的详细介绍可参阅文献[4-6]。由于太阳系行星历表误差难于建立模型，本文采用JPL公布的DE421，原子时采用BIPM发布的TT（2011）。

1 单引力波源引起的脉冲星计时残差

单引力波源分为两种：一种是单一持续的引力波源，比如双致密星的相互绕转；另一种是单一瞬时的引力波源，比如双致密星的并合。对于脉冲星计时阵所响应的引力波频率10-9～10-7Hz，致密星须为超大质量黑洞，特征质量一般为109M⊙，其中M⊙为太阳质量。下面分别讨论这两种引力波。

1.1单一持续引力波源

相互绕转的超大质量双黑洞在并合前的105～106a，其演化效应可忽略，由于辐射周期性的引力波离地球非常遥远，例如blazar OJ287，离我们有Gpc之遥，而脉冲星距离地球约Kpc的量级，因此引力波完全可视为平面波。于是，引力波在横向无迹规范下可写为[7]

式（11）和（12）中的αp和δp分别是脉冲星的赤经和赤纬。下面我们利用Tempo2模拟一个单一持续引力波源，其位置为：α=60°,δ=40°，频率为ν=10-8Hz，幅度为h0=10-13，这里需要说明的是，根据脉冲星计时阵对单引力波的限制，h0不应大于1.7× 10-14[11]，但是为了在计时残差中获得更加清晰的引力波信号，我们的取值提高了一个量级。另外，我们选取噪声较小的3颗脉冲星：J0437-4715，J1713+0747和J1909-3744。脉冲星的模型参数取自Parkes望远镜释放的数据“TheParkesPulsarTimingArray（PPTA）DataRelease1”。我们先分别模拟3颗脉冲星的TOA，测量误差引起的白噪声均方根设定为100 ns，然后加入模拟的引力波作为全局参数拟合。图1左列是3颗脉冲星加入引力波后拟合前的计时残差，这里计时残差的单位μs代表微秒。我们可以清楚地看到正弦波信号，而且不同脉冲星引力波引起的计时残差不同，这是因为引力波引起的计时残差不但与引力波的幅度有关，还与引力波与脉冲星的空间位置分布有关，这一点可以从式（9）中看出。图1右列是拟合后的残差，基本是白噪声，因此如果有如此强的单一持续引力波，多颗毫秒脉冲星的全局拟合很容易检测到它。需要指出的是，我们在上面的模拟中没有考虑“脉冲星项”，因为Tempo2对引力波信号作全局拟合时没有考虑“脉冲星项”，原因是引力波信号在计时残差中体现的相关性主要在“地球项”。忽略“脉冲星项”不会对结果造成很大影响。为检验Tempo2模拟单一持续引力波源的可靠性，我们将模拟结果和根据式（9）解析计算的结果做了比较。当然，作解析计算时我们也忽略掉“脉冲星项”。图2中的虚线是解析计算的结果，数据点及误差棒是模拟的结果，而实线是模拟数据作的曲线拟合。从图2中可以看出，解析结果和模拟结果在周期上符合得很好，但是在幅度上大约有26%的误差。

图1　单一持续引力波源引起的计时残差

图2　J0437-4715的模拟结果与解析计算比较

1.2单一瞬时引力波源

超大质量双黑洞由于持续辐射引力波，使得双黑洞系统的能量减少，因而绕转周期变得越来越短，彼此相距越来越近，最终双黑洞发生并合事件，产生幅度巨大的爆发性引力波。由于这种爆发性引力波包含非振荡项，导致了时空的永久性变形，因此其波形呈现永久性和非振荡性，通常被称作引力波“记忆”（GWM）[12]。引力波“记忆”效应可瞬间改变脉冲星的视周期，如果被脉冲星计时阵同时监测的话，脉冲星视周期改变依据脉冲星空间分布呈四极性。需要指出的是，引力波“记忆”与脉冲星中常见的“glicth”现象在脉冲星计时数据中表现不同，主要有以下几点区别：1）“glitch”在改变脉冲星自转频率ν的同时，通常会伴随频率一阶导数ν˙的改变，而引力波“记忆”只改变ν；2）“glitch”在改ν变后，过一段时间会恢复到原来的水平，而引力波“记忆”是永久性改变；3）目前的观测中除了磁星中观测到一个“anti-glitch”之外，所有的“glitch”只会增大ν，而引力波“记忆”可能使ν增大也可能使ν减小；4）“glitch”大都发生在年轻的脉冲星中，在毫秒脉冲星中很少发生，但引力波“记忆”产生的ν改变在毫秒脉冲星和正常脉冲星中都会出现；5）如果是脉冲星计时阵观测，则引力波“记忆”使得所有脉冲星的ν都发生改变，而“glitch”属于脉冲星的独立事件，只能在个别脉冲星中观测到。

下面我们讨论一下引力波“记忆”如何影响脉冲星计时残差。Favata[13]指出，对于圆轨道而言，引力波“记忆”只有“+”极化的分量，因此引力波“记忆”可以用一个阶跃函数表示[14]：

式（15）中，ε=ΔM/M 为贡献引力波“记忆”效应的质量比例，约为0.07[13]，μ=m1m2/（m1+m2）为双黑洞的约化质量，Dc为双黑洞的共动距离。则由引力波“记忆”引起的计时残差仍然可以用式（8）计算。下面我们利用Tempo2模拟。首先设定一些参数：hm=10-13，α=60°，δ=40°。同样选取上面的3颗毫秒脉冲星。图3给出了自转参数拟合前（左列）和拟合后（右列）的计时残差。很明显，3颗星拟合前的计时残差在同一时间有线性变化，表示改变了视周期，但是改变的方式及大小都与脉冲星的位置分布有关。

图3　引力波“记忆”引起的计时残差

2 随机引力波背景引起的脉冲星计时残差

随机引力波背景的成因主要分为3种：宇宙中多个双星系统的绕转和并合形成的随机背景、暴涨后期形成的残余引力波背景以及宇宙弦模型预言的引力波背景。这3种背景引力波的特征强度谱在脉冲星计时阵的响应频率范围内，都可以简单得表示为幂率形式：

式（18）中，σ2为由随机引力波背景引起的脉冲星计时残差的方差。目前，Shannon等人利用PPTA的观测数据将A限制为：A＜2.7× 10-15[17]。由式（16）～（18）我们可以得出，P（ ν）∝ν2α-3，σ2∝ν2α-2。下面我们只讨论两种背景引力波，即多个超大质量双黑洞系统产生的随机引力波和暴涨产生的残余引力波。为方便叙述，我们分别称作I型引力波背景和II型引力波背景。对于I型背景，谱指数为α=-2/3；而II型背景，谱指数为α≈-1。下面我们以计时红噪声较小的脉冲星J1713+0747为例，在其观测数据上加上模拟的背景引力波，分析它的计时残差变化。根据PPTA释放的数据，利用Tempo2可以拟合出J1713+0747的视差，得到其视差距离为1.12Kpc。在此给出两种情况的拟合后方均根残差。为更加明显地展示引力波背景引起的计时残差，我们把引力波的幅度比上限值提高一个量级，A=2× 10-14。引力波源数设定为1000，频率范围取为：下限νlow=0.01/Tspan，上限νlow=1/（1d），其中Tspan为数据总的时间跨度，1d代表一天。引力波谱指数分别取为α=-2/3和α=-1。如图4所示，左列是拟合前的计时残差，右列是拟合后的计时残差；第1行为没有引力波的情况，第2行为I型背景引力波，第3行为II型背景引力波。我们可以看到，在模拟参数和条件相同的情况下，不同引力波背景产生的计时残差不同，与其谱指数有关。II型引力波背景比I型引力波背景产生的计时残差要显著许多。因此，在相同强度下和相同源数下，II型引力波背景更有希望被探测到。当然，由于，ng为引力波源数，引力波引起的计时残差的均方根会随ng的增大而降低。从拟合前的计时残差可以看出，引力波背景对脉冲星的视周期影响非常明显，因为计时残差基本呈线性变化趋势。拟合后的计时残差呈现出类似低频红噪声的特征。但多颗脉冲星的计时残差会呈现出Hellings-Downs曲线的相关性。在这里我们不做讨论了，详细原理可参见文献[18]。

图4　引力波背景引起的脉冲星J1713+0747的计时残差

3 结语

本文利用脉冲星计时处理软件Tempo2模拟了3种形式的引力波对脉冲星计时残差的影响。单一持续引力波会在脉冲星计时噪声中留下周期性正弦波信号（对圆轨道而言），单一瞬时引力波可瞬间改变脉冲星的视周期，而引力波背景会在脉冲星拟合后的计时残差中呈现低频红噪声的特征。由于不同形式的引力波在脉冲星计时残差中表现出各自的特征，因此比较容易从计时残差中分离出来，从而得到引力波源的物理特征和属性。另一方面，从脉冲星时间标准或者脉冲星深空自主导航来讲，引力波却是我们数据处理中的噪声，这类似于脉冲星本身的计时红噪声，当然，计时红噪声无相关性，但引力波有相关性。更加清楚地认识并估计这些噪声，有利于建立综合脉冲星时间尺度，也有利于提高深空脉冲星导航的精度。因此研究计时噪声本身是有重要意义的。

[1]JENET F A,HOBBS G B,LEE K J,et al.Detecting the stochastic gravitational waves background using pulsar timing[J].ApJ., 2005,625:123-126.

[2]CHAMPION D J,HOBBS G B,MANCHESTER R N,et al.Measuring the mass of solar systemplanets using pulsar timing[J].ApJL,2010,720:201-205.

[3]HOBBS G,COLES W,MANCHESTER R N,etal.Developmentofpulsarbasedtimescale[J].MNRAS,2012,427:2780.

[4]HOBBS G,EDWARDS R T,MANCHESTER R N.TEMPO2:anewpulsar-timingpackage-I.Anoverview[J].MNRAS,2006,369: 655-672.

[5]EDWARDS R T,HOBBS G B,MANCHESTER R N.TEMPO2:a newpulsar-timing packageII.Thetiming model and precision estimates[J].MNRAS,2006,372:1549-1574.

[6]HOBBS G,JENET F,LEE K J,etal.TEMPO2:anewpulsar timingpackageIII.Gravitational wavesimulation[J].MRNAS,2009, 394:1945-1955.

[7]LEE K J,WEX N,KRAMER M,et al.Gravitational wave astronomy of single sources with a pulsar timing array[J].MNRAS, 2011,414:3251-3264.

[8]TONG Ming-lei,ZHAO Cheng-shi,YAN Bao-rong,et al.Pulsar timing residuals due to individual non-evolving gravitational wavesources[J].Res.Astron.Astrophys.,2014,14:390-400.

[9]WAHLQUIST H.The Doppler response to gravitational waves from a binary star source[J].Gen.Rel.Grav.,1987,19: 1101-1113.

[10]ANHOLMM,BALLMER S,CREIGHTON J D E,et al.Optimal strategies for gravitational waves stochastic backgroundsearches inpulsartimingdata[J].Phys.Rev.D,2009,79:084030.

[11]ZHU Xing-jiang,HOBBS G,WEN Lin-qing,et al.An all-sky search for continuous gravitational waves in the parkes pulsar timingarraydataset[J].MNRAS,2014,444:3709-3720.

[12]MADISON D R,CORDES J M,CHATTERJEE S.Assessing pulsar timing array sensitivity to gravitational wave bursts with memory[J].ApJ,2014,788:141.

[13]FAVATA M.Nonlineargravitational-wavememoryfrombinaryblack holemergers[J].2009,ApJ,696:L159.

[14]WANG J B,HOBBS G,COLES W,et al.Searching for gravitational wavememory bursts withtheParkes Pulsar Timing Array[J].2015,446:1657.

[15]CORDES J M,JENET F A.Detectinggravitational wavememorywithpulsartiming[J].ApJ,2012:752:54.

[16]JENET F A,HOBBS G B,et al.Upper bounds onthelow-frequency stochastic gravitational wavebackground frompulsar timing observations:currentlimitsandfutureprospects[J].ApJ,2006,653:1571-1576.

[17]SHANNON R M,OSLOWSKI S,DAI S,et al.Gravitational-wave limits frompulsar timing constrain supermassive black hole evolution[J].Science,2013:343-334.

[18]HELLINGS R W,DOWNS G S.Upper limits on the isotropic gravitational radiation background frompulsar timing analysis[J].1983,265:39-42.

Simulations and analyses of pulsar timing residuals induced by gravitational waves

TONG Ming-lei1，2，3，DING Yong-heng1，4，ZHAO Cheng-shi1，2，GAO Feng1，4，GAO Yu-ping1，2
（1.National Time Service Center，ChineseAcademy of Sciences，Xi′an 710600，China；2.Key Laboratory of Time and Frequency Primary Standard，National Time Service Center，ChineseAcademy of Sciences，Xi′an 710600，China；3.Key Laboratory for Research in Galaxies and Cosmology，ChineseAcademy of Sciences，Hefei 230026，China；4.University of ChineseAcademy of Sciences，Beijing 100049，China）

In order to separate the pulsar timing residuals induced by gravitational waves（GWs）from the total pulsar timing residuals，we introduced three kinds of GWs（single continuous GWs，single instantaneous GWs and stochastic gravitational wave background）and their influences on pulsar timing residuals.Moreover，we did some simulations using the processing software Tempo2 for pulsar timing，and gave some analyses and comparisons for the pulsar timing residuals obtained from the three kinds of the simulated gravitational waves.This work is of some significance in the construction of pulsar time scale and their applications.

pulsar；gravitational wave；timing residual

TN966+.7

A

1674-0637（2015）01-0044-08

10.13875/j.issn.1674-0637.2015-01-0044-08

2014-12-09

国家自然科学基金资助项目（11103024，11373028，11403030）；中国科学院“西部之光”人才培养计划西部博士资助项目（2011BS10）；中国科学院“西部之光”人才培养计划西部联合学者资助项目（2010LH02）；中国科学院星系宇宙学重点实验室开放课题资助项目（14010205）

童明雷，男，博士，副研究员，主要从事脉冲星计时与应用研究。