亲手复算HERA上限，并经窗函数比较模型

你需要知道什么是带功率、窗函数，以及规则 $U=\max(\hat\Delta^2,0)+2\sigma$。课程章节把时延功率变成带功率讲解了这些概念。

文件及其性质

所有文件放在同一个文件夹里，都可以单独下载。

摘录取自HERA合作组公开发布的HERA_2021_pspec_data.zip中的三个UVPSpec文件，下载于2026年11月12日，SHA-256为1c98d163...388d2c。论文写明，其表3、表4和图17、图18背后的数据是公开的。HERA Collaboration 2022, Data section 发布页没有写许可协议，所以这份摘录仅供教学和核验，并请使用者引用原论文。用h5py和numpy读取HDF5文件的提取脚本保存在本站源代码仓库里，其中从时延到 $k$ 的换算按文件里记录的git版本，照hera_pspec原样写出。hera_pspec conversions.py

模型文件不是HERA的产品，也没有拟合HERA数据。它就是实验室里的线性预测 $\Delta^2_{21}=\bar T_{\mathrm b}^2(1+\mu^2)^2\Delta^2_m$，取 $\mu=1$，用Eisenstein-Hu无振荡转移函数和Planck 2018参数。Eisenstein and Hu 1998 加热历史在 $z=7.93$ 处 $\bar T_{\mathrm b}=16.2\,\mathrm{mK}$，冷历史为 $-221\,\mathrm{mK}$。

运行

在放文件的文件夹里执行：

python3 -I recompute_window_comparison.py hera-h1c-idr2-window-extract.json model-spectra.json my-output.json

这个参数让Python不从当前文件夹导入任何东西。脚本会打印：

redshift 10.38
  k=0.256  U=(95.74)^2 mK^2  heated/U=0.0032  cold/U=0.375
  k=0.384  U=(178.23)^2 mK^2  heated/U=0.0014  cold/U=0.163
  k=0.512  U=(296.30)^2 mK^2  heated/U=0.0006  cold/U=0.076
redshift 7.93
  k=0.192  U=(30.76)^2 mK^2  heated/U=0.0168  cold/U=3.132
  k=0.320  U=(49.60)^2 mK^2  heated/U=0.0110  cold/U=2.063
  k=0.448  U=(113.39)^2 mK^2  heated/U=0.0029  cold/U=0.542
  k=0.576  U=(148.85)^2 mK^2  heated/U=0.0021  cold/U=0.389

然后把my-output.json和expected-output.json对照，只有Python版本一项允许不同。参考运行用的是Python 3.9.6，模型文件由Node v26.11.1导出。

如果某一行窗函数之和与一相差超过 $10^{-6}$，或模型网格与窗函数网格不一致，脚本会报错停下。

第一步：重建上限

文件里的带功率排在 $\Delta k=0.064\,h\,\mathrm{Mpc^{-1}}$ 的网格上。相邻的箱彼此相关，所以HERA的公开笔记本建议：要么每个频带只取最低的那个上限，要么隔一个箱取一个。脚本照笔记本的做法隔箱取值，第一频带从 $0.256$ 开始，第二频带从 $0.192$ 开始。HERA Public Data Release 1, notebook

功率单位为 $\mathrm{mK^2}$，$\sqrt U$ 单位为 $\mathrm{mK}$。两个标志性上限，$z=10.4$ 的 $(95.74)^2$ 和 $z=7.9$ 的 $(30.76)^2$，正是论文印出的数值。HERA Collaboration 2022, §4.4 精确复现它们，同时检验了摘录和规则。

第二步：让模型穿过窗函数

对每个保留的箱，脚本把HERA完整窗函数网格上的模型值乘以这个箱的窗函数行，再求和。它也给出箱中心处的模型值，方便你看窗函数到底起多大作用。在 $z=7.9$，冷历史在前两个箱超过上限，分别是 $3.13$ 倍和 $2.06$ 倍；加热历史在每个箱都不到上限的 $0.017$。在 $z=10.4$，两种历史都没有被排除。

这能说明什么、不能说明什么：完全耦合、从未加热的气体，即便只算密度起伏，在 $z=7.9$ 也被HERA第一天区的数据排除了。至于自旋温度或电离起伏很强的模型，这个模型没有包含这些项，所以对它们什么也说明不了。

五种会改变结论的错误

脚本有意算出每一种错误，并把结果写在正确值旁边。

1. 保留负的估计值。$k=0.320$ 处估计值为 $-516.0$，直接加 $2\sigma$ 得到 $1943.7$，而不是 $2459.7\,\mathrm{mK^2}$，上限比数据所能支持的严了21%。
2. 只加 $1\sigma$。$k=0.192$ 处上限从 $946.4$ 降到 $579.0\,\mathrm{mK^2}$，这就不再是95%的上限了。
3. 不用窗函数。$k=0.192$ 处冷模型在箱中心为 $2893$，经窗函数为 $2964\,\mathrm{mK^2}$；加热模型为 $15.49$ 对 $15.87$。这里差别不大，但它让每一次比较都朝同一个方向偏。
4. 把 $h$ 换算了两次。HERA的 $k$ 已经以 $h\,\mathrm{Mpc^{-1}}$ 为单位，再除一次 $h=0.6736$，每个点都挪到 $1.48$ 倍大的 $k$ 上，$0.192$ 处经窗函数的冷模型从 $2964$ 升到 $4604\,\mathrm{mK^2}$，高了55%。
5. 把每个原生箱都当成独立的。摘录里带有相邻箱之间的协方差，相关系数在第二频带为 $0.25$ 到 $0.31$，在第一频带为 $0.32$ 到 $0.33$。把七个原生箱当作独立数据合并，会高估数据的约束力。

单位、约定与这次核验的边界

文件里的 $k$ 采用HERA的宇宙学参数，$H_0=67.27\,\mathrm{km\,s^{-1}\,Mpc^{-1}}$。HERA Collaboration 2022, §4 模型采用Planck 2018，$h=0.6736$。两者 $h$ 相差0.13%，对模型点位置的影响远小于 $0.031\,h\,\mathrm{Mpc^{-1}}$ 的窗函数宽度，这里忽略不计，并记录在清单文件里。

估计值是交叉功率谱的实部，单位 $\mathrm{mK^2}$。误差是文件中协方差对角元的平方根；既然七个上限都与论文一致，它们就是论文表3和表4背后的误差棒。摘录里也有第二、第三天区，但比较只用第一天区，和论文的标志性上限一致。这不是似然分析。严格的约束要把相关的上限放进统计模型里处理。

检查清单

• 运行前，对照MANIFEST.json核对每个文件的SHA-256。
• 用python3 -I运行，并把输出与expected-output.json对照。
• 对照论文确认七个上限，尤其是两个标志性数值。
• 负的估计值置零，加 $2\sigma$，隔箱取值。
• 让模型穿过窗函数行，$k$ 只按 $h\,\mathrm{Mpc^{-1}}$ 换算一次。
• 画图时，HERA的数值标为真实观测，模型数值标为理想参考。