复算点透镜放大因子，再检验一次透镜搜索

先读点质量透镜放大因子的精确解和匹配滤波那一章。这里不用任何探测器数据，信号和噪声全是合成的。

下载文件

把下面八个文件放进同一个空文件夹：

• README.md：运行说明和失败案例
• point_lens_reference.py：闭式解和两像极限
• independent_fresnel_integral.py：计算 $F$ 的第二条路线
• signal_injection.py：噪声、啁啾、注入和搜索
• reproduce.py：全部重跑并比对
• injection-inputs.json：每次注入的设置和随机种子
• expected-outputs.json：本页引用的全部数字
• provenance.json：来源、常数和 SHA-256 校验和

然后运行：

python3 -m pip install numpy scipy mpmath
python3 reproduce.py

脚本会打印与保存结果的最大差异。所有放大因子相差不超过 $2\times10^{-6}$、所有注入中位数相差不超过 $10^{-6}$ 时，退出状态为 0。我们的环境是 Python 3.9.6、numpy 2.0.2、scipy 1.13.1 和 mpmath 1.3.0。换了版本，随机噪声的末几位可能变化，严格比对就会失败；这时用 python3 reproduce.py --write 保存你自己的结果，中位数应当和我们的接近。

单位与约定

太阳质量参数取 IAU 2015 名义值 $1.3271244\times10^{20}\,\mathrm{m^3\,s^{-2}}$，再除以 $c^3$。

第一步：闭式解

point_lens_reference.py 用 mpmath 以 60 位有效数字计算 Takahashi 和 Nakamura（2003）的式 (17)，以及他们的两像公式式 (18)。式 (17) 完全按印出的形式使用。Takahashi 和 Nakamura 2003，第 2 节

第二步：完全不用闭式解的第二条路线

independent_fresnel_integral.py 从透镜平面积分重新出发，把透镜平面上的点按到达时间分组。设 $A(\tau)$ 是 $T<\tau$ 的区域面积，则

第二个等号用到 $\int_0^\infty e^{iw\tau}\mathrm d\tau=i/w$。$I(\tau)$ 是透镜对一个尖脉冲的响应，没有透镜时等于 1，也就是 2003 年论文式 (16) 的时域图像。

脚本分四步算 $I(\tau)$：

1. 对从透镜出发的每个方向，在 $\ln r$ 上用二分法找到 $T=\tau$ 的位置，射线上极小值两侧各找一个。
2. 用 240 点高斯–勒让德求积对方向积分，先做变量替换，消掉等值线闭合处的平方根奇点。
3. 到达时间超过鞍点后，把求积节点集中到透镜另一侧，那里等值线收窄。
4. 对列表中的 $I(\tau)$ 按分段线性函数做精确傅里叶积分，$\tau=400$ 以外的尾部用 $c_1/\tau+c_2/\tau^2$ 拟合补上。

这条路线完全不碰伽马函数和超几何函数，所以实部和虚部都对上，检验的就不只是式 (17) 的模，还有它的相位约定。

结果。 在 $y\in\{0.1,0.3,1,3\}$、$w\in\{0.1,1,3,10,30\}$ 上，最大相对差是 $4.2\times10^{-4}$，出现在 $y=0.3$、$w=30$。每个 $y$ 拟合出的尾部系数都是 $c_1=0.499$。

全部 20 行都在 expected-outputs.json 里。

到达时间里的常数在这里很要紧。如果照 2003 年论文印出的式 (12) 在 $T$ 里写 $+\phi_m(y)$，菲涅耳路线得到的 $F$ 会多出相位 $e^{2iw\phi_m}$，就和式 (17) 对不上了。论文正文说最早到达时间取为零，这要求写 $-\phi_m$。我们按这个做，除此之外不对论文下任何结论。

第三步：两个极限

这些检验用的是闭式解，它已经在第二步得到验证。

• 低频。 $\lvert F-1\rvert$ 应当几乎与 $w$ 成正比地减小，因为有 $w\ln w$ 项，减得稍慢一点。$y=1$ 时，$w=10^{-2}$、$10^{-3}$、$10^{-4}$ 对应 $2.2\times10^{-2}$、$3.3\times10^{-3}$、$4.4\times10^{-4}$。$y=0.1$、$w=10^{-4}$ 时是 $5.13\times10^{-4}$，一阶公式 $\pi w/4+(iw/2)[\ln(w/2)+\gamma-2\phi_m]$ 给出同样的值。
• 高频。 与两像公式的差除以 $\lvert F\rvert$，应当随 $w$ 增大而减小。$y=1$ 时从 $w=30$ 的 $2.6\times10^{-3}$ 降到 $w=1000$ 的 $2.6\times10^{-5}$。$y=0.3$ 时，$w=30$、100、300、1000 依次是 $3.1\times10^{-2}$、$6.5\times10^{-3}$、$8.7\times10^{-3}$、$1.1\times10^{-3}$。100 到 300 之间的回升是真实的：剩余的差本身随 $w$ 振荡。

第四步：注入与搜索

signal_injection.py 生成一段 8 s、每秒 1024 个采样的数据：

1. 噪声。 按 Ajith（2011）式 (4.7) 着色的高斯噪声，每个频率格的实部和虚部标准差都是 $\sqrt{S_n/(4\Delta f)}$。Ajith 2011，第 IV 节
2. 信号。 $10+10\,M_\odot$ 的最低阶啁啾，从 20 Hz 到 219.9 Hz，缩放到无透镜信噪比 20；透镜注入再乘上精确的 $F$。Cutler 和 Flanagan 1994，式 41–42
3. 无透镜搜索。 黑洞质量从 9.6 到 $10.4\,M_\odot$ 取 17 个值，每个都对幅度、到达时间和相位取最大。
4. 透镜搜索。 透镜质量在 10 到 $1500\,M_\odot$ 之间取 36 个值，$y$ 在 0.1 到 3 之间取 9 个值；然后对每个 $y$，在最佳质量的两个相邻格点之间再细扫 48 个质量；最后对最好的两组做有界的一维优化。透镜模板保持真实的黑洞质量。
5. 统计量。 $\Delta\ln L$，即最佳透镜模板与最佳无透镜模板的对数似然比之差，报告门槛取 4。

加噪声之前，脚本先检查搜索能否找回注入的透镜。两种透镜情形都找回了：透镜质量误差在 0.001% 以内，$\Delta\ln L$ 和精确模板的值相差不到 $10^{-8}$。

结果。

清楚的透镜不会认错。弱透镜像抛硬币，而纯噪声的尾巴能伸到弱透镜也会给出的数值。

有意保留的失败案例

• 格点太粗。 第一版只在 36 个质量格点上搜索，相邻格点相差约 15%。$y=2$ 时延迟是 90 ms，相邻质量之间 200 Hz 处的条纹相位会差好几个弧度，所以即使没有噪声，最好的格点模板也只有 $\Delta\ln L=0.73$，而不是 6.03。那一版弱透镜的中位数 2.79 基本是噪声。第四步的细化修正了这个问题，无噪声检查也能再次发现它。
• 放大因子表太粗。 第一版 $F$ 表在 $w$ 上只有 1500 个点。$y=3$、$w$ 约为 41 时，每个条纹只采到约四个点，插值误差达到 3.8%。改成 6000 个点后降到 0.24%。
• 假透镜。 120 次纯噪声试验中有 14 次越过门槛。最大的一次是种子 3058，数值 7.97，偏好 $M_{Lz}\approx1150$、$y=1.5$。
• 没有包含的东西。 没有自旋、不等质量、并合、铃宕、探测器天线方向图、突发干扰，也没有透镜周围的星系；$y$ 只在格点上搜索。这些数字用来说明方法，不是任何真实搜索的灵敏度。

相信一次重跑之前，逐项核对

• python3 reproduce.py 退出状态为 0，或者你自己用 --write 得到的结果只在与噪声有关的末位上和我们不同。
• provenance.json 里的 SHA-256 校验和与你的文件一致。
• 菲涅耳路线和式 (17) 在实部和虚部上都一致，而不只是模一致。
• 低频时 $\lvert F-1\rvert$ 大致按 $w\ln w$ 减小，高频时两像误差总体在减小。
• 读任何带噪声的结果之前，无噪声搜索先找回了注入的透镜。
• 每个透镜结果旁边都报告了纯噪声的分布，并且用同一个门槛。