用已发表的能流表复算一条 EMRI 频率演化

要复算的结果

三条圆赤道旋近， $q=10^{-5}$ ， $M=10^6\,M_\odot$ ：零自转从 $r=10$ 到 ISCO 外侧，自转0.9从 $r=10$ 到 $r=6$ ，自转 $-0.9$ 从 $r=20$ 到 $r=10$ 。每一条，脚本都给出用表中能流和用四极能流的时间与波周数，以及九个点的频率轨迹。

情形	周数（表中能流）	年数（表中能流）	周数（四极能流）	比值
$a=0$ ， $10\to6.003$	$112\,981$	$1.417$	$112\,045$	$0.9917$
$a=0.9$ ， $10\to6$	$260\,673$	$3.157$	$222\,719$	$0.8544$
$a=-0.9$ ， $20\to10$	$971\,669$	$31.47$	$970\,777$	$0.9991$

第一步：冻结的能流表

Fujita 和 Tagoshi 的表VII给出圆赤道轨道到 $\ell=6$ 为止的引力波光度，自转从 $-0.9$ 到 $0.9$ 共七种。输入文件收录 $r=6$ 到 $r=20$ 各行，是从 arXiv 第1版手工抄录的，表注原文一并保存。

表里没有可以直接引用的归一化说明，所以“按 $(\mu/M)^2$ 归一、到达无穷远的能流”这个解读是我们的推断。有两项检验支持它： $r=6$ 处逐个模式的表III加起来，与表VII的值相差不到 $5\times10^{-5}$ ；在 $r\ge50$ 的远处，表中数值与 Tagoshi–Sasaki 后牛顿无穷远能流相差不到 $2\times10^{-4}$ 。

对应阅读：Fujita & Tagoshi (2004) · Tagoshi & Sasaki (1994)

第二步：插值一个平缓的修正因子，而不是能流本身

整张表里能流跨了三个数量级，所以脚本插值的是表中能流与四极能流之比。这个比值始终在0.8到2之间，随半径平缓变化。脚本用最近的四个表列半径做三次插值，6到20以外的半径一律拒绝。

能量斜率 $E'(r)$ 由 Kerr 圆轨道精确能量做四阶有限差分得到。两个积分都在 $\ln r$ 上用辛普森公式，分2000段。

t=\int_{\ln r_1}^{\ln r_0}\frac{r\,E'(r)}{\mathcal F(r)}\,d\ln r,\qquad N=\int_{\ln r_1}^{\ln r_0}\frac{\Omega_\phi(r)}{\pi}\,\frac{r\,E'(r)}{\mathcal F(r)}\,d\ln r,\qquad \mathcal F=\tfrac{32}{5}\,q\,\Omega_\phi^{10/3}\times\frac{\mathcal F_{\rm table}}{\mathcal F_{\rm quad}}

在

\ln r

上积分，能把取样点均匀地放在被积函数变化最快的地方。

第三步：运行并比较

下载文件夹，先用校验和确认文件完整，用输入文件和一个新的输出文件名运行脚本，再对照你的输出和预期输出。在生成预期文件的机器上（macOS，Python 3.9.6），重新运行逐字节一致。换一台机器，预计相对差在万亿分之一左右。

时间用 $GM_\odot/c^3=4.925490947\times10^{-6}\,\mathrm s$ 换成秒，用儒略年31,557,600秒换成年。轨迹里的波频率是 $f=2\Omega_\phi/(2\pi)$ ，单位赫兹，对应 $M=10^6\,M_\odot$ 。

第四步：故意改坏

输出里记了三处故意的改动，以及每处让零自转周数变了多少。修正因子改用线性插值：−0.54%。辛普森公式从2000段减到4段：+0.088%。把 Kerr 能量换成牛顿束缚能：+99%。

最后一个数最值得记住：靠近 ISCO 时，束缚能和能流一样要紧。要公平地检验一个能流模型，轨道本身必须是精确的。

结果说明不了的事

拒绝外推。自转0.9的 ISCO 在 $2.32$ ，远低于表的第一行，所以这一情形停在 $r=6$ 。输出里把 $r=4$ 的请求列为已拒绝。

缺少的多极。 $r=6$ 处零自转能流里， $\ell=6$ 模式占0.29%，略去的 $\ell\ge7$ 模式那里大约占0.1%。这是顺着趋势的外推，不是计算。

缺少的物理。视界能流、偏心和倾斜轨道、共振，以及所有后绝热项都没有算进去。后绝热相位修正约一弧度，所以这些周数只在质量比的领头阶上是对的。

对应阅读：Fujita & Tagoshi (2004) · Barack & Pound (2019)

下载并运行

把六个文件放进同一个文件夹，运行下面两条命令。第二条会在预期输出旁边写出你自己的结果。

shasum -a 256 -c SHA256SUMS python3 recompute_frequency_evolution.py inputs.json my-output.json

recompute_frequency_evolution.py · 脚本，只用 Python 标准库
inputs.json · 情形、单位与冻结的能流表
expected-output.json · 预期输出，Python 3.9.6
README.txt · 单位、模型与失败情形
SHA256SUMS · 每个文件的 SHA-256 校验和
LICENSE.txt · 使用许可与署名

在实验室里试一试

与实验室的四极周数对照

在实验室里设自转为零、偏心率为零、 $p=10$ 。
读出“剩余波周数”。实验室保留三位有效数字，与上表第一行四极列的112,045比较。

实验室从你设定的轨道大小开始计数，只用四极能流。

练习与答案

哪个情形里表中能流与四极能流差得最多？为什么是它？

自转0.9、从 $r=10$ 到 $r=6$ ，比值0.854。在这段范围里四极能流处处偏大15%到19%，误差不变号，也就没法抵消。零自转时误差在 $r=8$ 附近变号。

把输入里零自转的终止半径改成5.9再运行，会发生什么？

脚本报错停下，指出这个半径并拒绝外推，因为 $r=5.9$ 低于表的第一行。即使表更大，它也落在零自转 ISCO $r=6$ 以内，那里没有稳定的圆轨道。

一手来源，以及每篇读了哪些部分

Barack & Pound (2019), Self-force and radiation reaction in general relativity, Rep. Prog. Phys., doi:10.1088/1361-6633/aae552
arXiv 第2版，2018年11月19日。第5.2节与式(72)至(74)：双时间尺度相位展开，以及绝热阶与后绝热阶各需要什么；第5.3节瞬时共振；第6.1节与式(77)至(78)，按轨道平均的平衡律，包括为负的视界能流。
Fujita & Tagoshi (2004), New numerical methods to evaluate homogeneous solutions of the Teukolsky equation, Prog. Theor. Phys. 112, 415
arXiv 第1版。表VII，七种自转下圆赤道轨道到 l = 6 的引力波光度；表III，r = 6M 处逐个模式的光度。单位归一化以及“到达无穷远的能流”这一解读是我们的推断，已用四极与后牛顿极限检验。
Tagoshi & Sasaki (1994), Post-Newtonian expansion of gravitational waves from a particle in circular orbit around a Schwarzschild black hole, Prog. Theor. Phys. 92, 745
arXiv 第1版。式(43)：到四阶后牛顿的无穷远能流与四极值之比；只用作能流表弱场一端的独立检验。

继续阅读

这一页直接连着哪些内容

当前页面用已发表的能流表复算一条 EMRI 频率演化

实验

本文用到小黑洞绕自转黑洞的极端质量比旋近

文章

极端质量比旋近 · 阅读路径