双黑洞系统的不同起源及其参数分布

2019-04-10吕振伟李木子郭越凡

天文学进展 2019年1期

吕振伟，李木子，郭越凡，张帆

(北京师范大学天文系，北京100875)

1 引言

物体的运动导致空间中质量能量分布的改变，从而引起时空曲率的微小变化；而时空的曲率又决定着物体如何运动。时空曲率的扰动如涟漪一般在宇宙中以光速传播且携带能量，这就是基于1916年爱因斯坦广义相对论所预言的引力波辐射。2016年2月，LIGO(Laser Interferometer Gravitational Wave Observatory)科学合作组织宣布位于美国Hanford和Livingston 的两台激光干涉引力波探测器首次成功捕捉到来自双黑洞并合的引力波信号GW150914[1]。这是引力波天文学史上具有里程碑意义的事件，为百年来物理学家们对于引力波是否为可测量的物理波的争论画上了句号，且再次验证了广义相对论的正确性。

目前高新LIGO[2]和高新VIRGO(Virgo interferometer)[3]已经升级完毕并投入到联合探测中，截止到2017年底(O1 和O2 已运行，2019年将开始O3 运行)，已有5 个双黑洞系统并合引力波信号[1,4–7]及1 个双中子星并合信号GW1710817[8]被HLV (Hanford-Livingston-Virgo)引力波探测器网络成功探测到。对双黑洞系统的质量估计在表1 中列出。为降低低频上的地震噪声，日本3 km 臂长的KAGRA 引力波探测器建于地下，且应用了冷冻镜面的技术，以改善热噪声的干扰，于2016年3月开机并进行了为期1 个月的试运行，预计2018年将正式投入观测[9]。2024年，LIGO-india 也会加入到全球地基引力波探测器网络中[10]。空间引力波探测项目如激光干涉空间天线(Laser Interferometer Space Antenna, LISA[11])，DECIGO[12]，太极(TAIJI[13,14])和天琴(TIANQIN[15])的优势是可在无需考虑地球曲率的情况下增加臂长，且避免了地球上的地震噪声和人为噪声的干扰，使其在低频上的灵敏度大大提升。LISA 的设计臂长为2.5×106km，工作频段在10−4∼10−1Hz，其科学目标为探测星系中心超大质量双黑洞系统并合、极端质量比双黑洞并合等低频引力波信号。

表1 已观测到的双黑洞并合引力波事件(M⊙表示太阳质量)

作为地基引力波探测器和空间引力波探测器的重点目标，双黑洞系统附近的物质通常比较稀薄，因此其绕转-并合- 铃宕的演化过程在传统天文观测的电磁波段上不可见。而在该过程中，系统能量以引力波的形式向外辐射，且由于引力波与物质之间的相互作用比较微弱，引力波信号可从被极强引力效应支配的近密双星并合发生区域逃逸，并且携带黑洞系统的直接信息，如双黑洞的质量、自旋、轨道离心率等。当高新LIGO 和高新VIRGO 在其设计灵敏度上运行或者第三代引力波探测器建设完成时，我们将能探测到双黑洞系统各阶段引力波辐射的更多细节，而高探测率带来的大量信号样本既是进一步研究黑洞物理的宝贵资源，又是对引力波信号处理和波源参数估计方法的全新挑战。

目前关于双黑洞系统起源的研究还不深入，而我们特别关注恒星级双黑洞系统，因该质量范围的黑洞在宇宙中大量存在，且有间接的光学观测可相互验证。普遍认为恒星级双黑洞系统的起源主要有大质量双星演化起源和动力学起源[16–18]这两种情况。不同起源的系统具有不同的黑洞自旋和质量分布，大质量双星演化起源会使得双黑洞自旋角动量方向和公转轨道角动量方向具有高度同向性，而动力学起源则会导致自旋朝向的随机分布。前身星核坍缩后的爆发一般是各向异性的，这可能会对自旋-轨道角动量方向夹角有重要影响[19,20]，因此大质量双星演化是否能如理论预言的一样，使自旋-轨道角动量同向还需要验证。我们仍然假设自旋-轨道角动量是同向的，并且有小角度弥散。如何利用引力波信号携带的波源信息进行参数估计，自旋-轨道角动量方向的同向性是否能用来推断系统的起源，回答这些问题对进一步理解双黑洞系统背后的物理至关重要。

本文主要面向刚入门引力波研究领域的学者和学生，介绍如何通过分析引力波信号推断波源参数来限制双黑洞系统起源与演化模型的方法。我们假设黑洞形成的天体物理过程决定了双黑洞质量和自旋等参数分布的大致轮廓；而利用引力波信号推断的波源质量自旋信息和动力学方法测量的结果来验证不同起源双黑洞系统质量、自旋分布的理论预言。综合考虑理论模型和实验观测的结果，就可对于不同起源双黑洞系统给出一个参数化的模型来描述黑洞质量分布与自旋分布。随后，我们对所构建的模型进行采样并模拟一个数据集，其中一部分为具有自旋-轨道同向性的大质量双星演化起源，另一部分为动力学起源。将该数据集作为LALInference[21]的输入(该软件提供了引力波参数估计的标准方法，见4.1 节)，我们得到质量和自旋等参数的后验分布。进而计算贝叶斯系数(Bayes factors)和总质量差等对模拟系统起源进行推断，得到的结果可验证利用自旋分布作为判据区分黑洞系统起源的有效性，回答数据集中不同起源所占的相对比例，以及不同起源总质量分布的差异等问题。

本文的结构如下，第2 章中介绍已有的各种理论演化模型，总结它们的异同；第3 章中，我们结合理论模型和观测结果，重新构建黑洞的参数模型；在第4 章中，进行采样并得到模拟数据集，从而利用LALInference 软件包将理论预言不同起源的双黑洞系统质量、自旋等参数的分布情况转化成可以与引力波观测结果直接比对的后验分布，进而达成限定诸如“不同起源双黑洞系统所占的比例”这类基本天体物理问题的目的。第5 章中，我们对本方法的结果进行讨论，总结现阶段工作中还需提升的部分及对未来工作的展望。

2 恒星级黑洞的观测和理论模型

根据引力波探测的结果，人们发现了另外一种恒星级黑洞质量分布规律。如首次探测到的双黑洞质量分别为29M⊙和36M⊙，并合后黑洞质量更是高达62M⊙[1]，以前从未探测到类似结果。而且，多数探测到的双黑洞并合事件的黑洞质量在20M⊙∼40M⊙之间[1,4–7]，这为黑洞以及双黑洞系统存在多种起源的观点提供了观测支持。

恒星级黑洞的理论模型主要有两种。一种是综合了恒星模型和超新星流体动力学模拟的结果，由Belczynski 等人[23]在2012年提出的大质量恒星演化的物理模型，该模型解释了瞬态低质量X 射线双星中黑洞的质量观测结果，并且合理地解释在中子星和黑洞之间的质量缺失问题。但对于黑洞质量大于15M⊙的情况，无法进行描述。另一种是脉冲对不稳定型超新星(pulsational pair instability supernovae)理论[24,25]，该理论可以解释引力波探测的大质量黑洞现象，这主要是由于大量的物质在坍缩之前被喷出。

Belczynski 等人[23]提出的大质量恒星演化的物理过程，对超新星爆发内部工作机制给出了严格的限制。特别是核坍缩型超新星是在恒星坍缩后100∼200 ms 爆发的，这意味着爆发是由10∼20 ms 内生长起来的不稳定性驱动的。或者说，如果将来的观测发现了缺失的黑洞，那表明这种不稳定性需要大于200 ms 的生长时间，这样才能形成质量为3M⊙∼5M⊙的致密天体。

脉冲对不稳定性超新星理论[24,25]可以解释引力波观测中的大质量黑洞现象。由于大量的物质在坍缩之前被喷出，这种不稳定性发生在大质量恒星(80M⊙∼150M⊙)演化的后期，伴随着高温低密度区引起的正负电子对热聚集，会对状态方程有重要的影响，也只有质量在80M⊙∼150M⊙范围内的恒星有足够高的熵可以产生这种不稳定性。由此产生的黑洞质量小于40M⊙，并且大多数黑洞的质量约为34M⊙，这与引力波观测的结果一致。对于质量在150M⊙∼250M⊙之间的原恒星，由于对不稳定性(pair instability)而不能形成黑洞。

3 不同起源的双黑洞系统的质量分布建模

目前恒星级双黑洞的起源主要有两种：大质量双星演化起源和动力学起源[16–18,26]。大质量双星演化发生在星系场(galactic fields)中，双星通过大质量双星演化形成双黑洞，具有自旋和公转轨道高度同向的特点。动力学起源发生在球状星团(globular clusters)等恒星密度高的环境中，自旋和轨道一般是随机的。下面我们先构造双黑洞系统中质量较大的主黑洞(m1)的质量分布模型，然后应用幂律分布构造质量较小的次黑洞(m2)的质量分布[27,28]。由于光学观测上黑洞质量在5M⊙∼15M⊙范围内，存在2M⊙∼5M⊙范围内的质量缺失，并且大多数在(7.8±1.2)M⊙范围内[22]。结合Belczynski 的理论模拟结果，可以假设小质量黑洞(5M⊙∼20M⊙)服从幂律分布(power-law distribution)，并且M5M⊙。由引力波探测得到的黑洞质量则大得多，主要分布在20M⊙∼40M⊙之间[1,4–7]；同时，利用脉冲对不稳定性超新星理论得出：大质量恒星(80M⊙∼150M⊙)演化得到的黑洞质量小于40M⊙，并且大多数为34M⊙。因此，可以用高斯分布来描述大质量黑洞(20M⊙∼40M⊙)的分布，两部分叠加得到[24]：

其中，λ是两部分(主黑洞和次黑洞)所占比例，H(mmax−m1)是截止函数，mmax是幂律分布的截止质量，mpp是高斯分布的平均质量，方差为σpp，Λ={α,β,λ,mmin,mmax,mpp,σpp}，mmin是幂律分布的初始质量，p(m1|Λ)是在Λ的条件下m1的条件分布。例如，可以取ΛA={1.5,3,0.2,3,35,38,1.5}和ΛB={1.5,2,0.1,3,35,35,1}，分别作为大质量双星演化起源(模型A)和动力学起源(模型B)的质量分布，用来进行MCMC 采样产生样本，如图1所示。

图1 模型A 和模型B 的质量分布

对于来自大质量双星演化的双黑洞系统，我们假设黑洞诞生时爆发的各向异性对自旋倾角的影响很小，它们的自旋轨道倾角在[0◦,10◦]内均匀分布或者符合高斯分布；而对于来自动力学起源的双黑洞系统，其自旋和轨道朝向是随机分布的。最后用马尔可夫链蒙特卡罗(Markov chain Monte Carlo, MCMC)或者多级抽样(nested sampling)[21]的方法进行采样，得到一组引力波事件。把这些事件对应的参数作为样本输入LALInference，可以得到相应的后验概率。

4 贝叶斯推理及结果

贝叶斯推理的特点是每次探测到的信号都可以用以前的后验概率作为先验概率，然后计算新的后验概率，通过连续叠加得到最终的后验概率。这是一种被天文学家广泛应用且行之有效的检验方式，可以有效地减少计算量，节约大量时间和计算成本，这对于计算量非常大的引力波数据分析工作至关重要。在接下来的数年内，我们会得到越来越多的引力波事件，如何精确地推断致密双星的起源，不同起源所占的比例，以及质量自旋等参数的分布情况，是需要解决的三个主要问题。

4.1 LALInference模拟引力波模板并进行贝叶斯推理

该软件是由C 语言编写，并且包含了用Python 写的后期处理工具，是处理引力波探测数据的标准方法，可以用来计算功率谱和给定噪声曲线的稳定高斯噪声，还可以使用所有的引力波模板；在给定模板和模型参数的条件下，得到数据的似然函数，以及参数空间的随机采样等。下面对LALInference 的基本原理进行简单介绍。

本文主要采用其参数估计部分，当输入一个引力波事件的数据时，LALIference 会比较输入的探测信号和内置的理论模板，得到关于该事件参数(质量、自旋等)的似然函数(likelihood function)，也就可得到在探测到该事件的条件下参数的条件分布，即参数的后验分布：

其中，θ包括了所有的模板参数，且∫p(θ)dθ= 1，d是探测到的数据集合。p(θ)是基于目前所有的信息假设的先验分布。p(d|θ)是似然函数，可以由探测信号的统计性质给出。p(d)被称为统计证据，是似然函数和先验相乘后对所有参数的求和，在参数估计时，p(d)被用来归一化。当我们输入不同的模板参数，LALInference 内置模板函数可以计算相应的波形模板，然后通过比较波形模板和探测数据得到似然函数。

一般我们假设探测器的噪声是稳定的高斯噪声，因此，噪声的似然函数可以用简单的高斯分布来描述：

图2 双黑洞并合引力波事件GW151226 的噪声功率谱

当探测数据中有引力波信号时，假设信号是h，那么n=d −h，包含引力波信号的似然函数：

下面我们由LALInference 相关函数得到的不同质量黑洞系统的引力波信号波形，进一步了解相应的频率和绕转速度变化曲线。我们在计算时，设双黑洞质量相等，即m1=m2。图3 和图4 仅给出引力波的h+偏振。

图3 不同质量的双黑洞系统引力波波形

图4 不同质量的双黑洞系统引力波波形

双黑洞质量相等(m1=m2)时，在相互绕转到并合时期的频率和绕转速度随时间变化情况如图5 和图6 所示。可以看到，在并合前数毫秒，频率和绕转速度指数增加。

图5 并合前引力波频率随时间变化

图6 并合前双黑洞绕转速度随时间变化

对于先验概率分布，假设起源A 的质量在10M⊙∼80M⊙之间，为均匀分布；起源B的质量在10M⊙∼80M⊙之间，呈高斯分布，方差σ2= 400。由此可以进一步计算后验分布，结果如图7 所示。

图7 质量分布采样

对于引力波源GW150914，假设质量自旋等参数都是均匀的先验分布，通过5 000 个样本计算后验概率。可以得到质量估计如果用更多的样本，就能得到更加严格的限制。某些参数组合的散点分布如图8 所示。

图8 质量m1 和m2 后验分布以及赤经和赤纬后验分布的散点图

4.2 推断引力波波源

假设探测到一个新的事件，通过参数估计可以得到波源各个参数的分布(贝叶斯估计得到的都是概率描述，概率论的观点得到的参数是确定值)。该并合事件的起源及其概率可以通过计算贝叶斯系数来进行推断[24]。假设HA是大质量双星演化的模型假设，HB是动力学起源的模型假设，d是模拟的引力波波形数据，则在模型HA的条件下，得到事件d的概率：

由于样本是服从p(θ|d)分布的，所以上式中的积分可以转化为求和：

式中，θ是所有的模型参数，对这些参数积分后得到的统计证据p(d|HA)是与具体参数无关的量。p(θi)是假设的先验分布，p(θi|H)的后验分布皆可通过LALInference[21]软件得到。模型A 和模型B 下，我们进一步比较引力波源的统计证据，得到：

其中，表示模型假设A 和B 的统计证据之比。当1 时，表示引力波信号来自模型A 的可能性比来自B 的可能性大很多。一般选取= 8 作为模型判断的标准，当8 时，我们就认为信号来自模型A。

4.3 推断不同起源所占比例

假如我们有大量探测事件，如何得到每种起源所占的相对比例？以及大概需要多少事件？我们可以模拟产生引力波信号样本，其中给定比例的事件来自大质量双星演化，其自旋和轨道朝向高度一致；另一部分来自动力学起源。利用贝叶斯推理，来自大质量双星演化的样本所占的比例，可以被有效地区分出来。也就是说，可以用来检验引力波是否可以用来区分自旋和轨道同向与否。这对于理解致密双星演化动力学过程至关重要。

假设起源A 所占的比例为CI，B 是CII，CI+CII=1，则后验概率为：

其中，

为便于计算，等式两边同时取对数，得到：

其中，p(CI|d)是在数据d的条件下，得到起源A 在所有事件中所占的比例CI的分布。p(CI)是给定的先验分布，p(dk|CI)是在CI条件下探测到事件k的概率。随着事件数的增加，可以得到p(CI|d)随事件个数变化的规律，即可估计出确定每种起源所占的相对比例的事件数，这对于双黑洞系统的起源研究和双黑洞并合率的研究至关重要。

4.4 质量分布的差异

我们假设不同起源的双黑洞具有不同的质量分布，同时也具有不同的自旋分布，然后使用自旋分布信息来区分不同起源的黑洞，检验不同起源质量分布的差异，例如，对两种起源的总质量进行比较。换言之，我们采用自旋表征的样本来对不同起源的质量分布进行估算。由金斯半径估算可得，场星系里由大质量双星演化诞生的双黑洞质量比球状星团里的要小一些，因此可假设大质量双星演化的黑洞比动力学起源的黑洞总质量小。通过对分布π(m|Λ,H)∝p(m1|Λ,H)·p(m2|m1,Λ,H)·Vobs(m)进行采样，从而得到不同起源的质量分布[24]。q=MI−MII的分布，其中MI=mI1+mI2,MII=mII1+mII2，p(q)由下式表示：

其中，

以上采用的是贝叶斯等级推理(Bayesian hierarchical inference)的方法，数据模拟可以通过LALInference 软件来完成。

5 总结与展望

本文介绍了引力波信号的后期处理，主要是用LALInference 模拟引力波模型并进行贝叶斯推理；然后利用贝叶斯公式推断引力波起源、不同起源的比例以及不同起源质量分布的差异检验公式，对贝叶斯等级推理的公式进行系统的推导；最后从理论上讨论了贝叶斯方法的可行性。

由致密双星并合产生的引力波信号包含了丰富的波源信息，包括致密双星的质量、自旋、轨道倾角、离心率等。我们主要研究质量和自旋两个参数在推断引力波源和内禀参数分布上的应用。事实上，引力波的波形中包含的参数高达15 个，当我们考虑更多的参数时，很难找到多参数分布的简单描述，计算量也会大大增加，实际操作受到很大的影响。这时可以考虑采用机器学习的方法，通过对构建的多层神经网络进行训练，可能会给出更加精确的解。并且，深度神经网络的方法虽然在训练深度网络模型时会耗费大量的时间，但是对于已经训练好的神经网络模型，几乎可以在瞬间对上述问题进行判断。以上的问题是未来引力波天文学的主要问题之一，我们希望随着观测事件的增加，能提出更好的理论模型来解释这些问题。